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ABSTRACT 


This  report  reviews  work  carried  out  to  identify  promising  techniques  for 
determining  the  major  acoustic  mechanisms  governing  scattering  from  the 
subbottom  in  shallow  water.  Two  approaches  were  investigated:  (1 )  the 
application  of  T-matrix  scattering  theory  to  scattering  from  geological  features  in 
the  subbottom,  and  (2)  the  coupled  differential  equation  theory  of  gradient  driven 
coupling  of  shear  and  compressional  waves.  T-matrix  theory  is  well  developed 
and  has  most  of  the  characteristics  needed  for  application  to  modeling  scattering 
from  the  subbottom.  T-matrix  theory  is  a  full  wave  approach,  has  few 
fundamental  restrictions,  includes  both  compressional  and  shear  waves,  applies 
to  multiple  structured  scatterers  of  quite  general  shape,  and  can  be  formulated 
for  a  waveguide  geometry.  The  basic  theory  of  gradient  driven  coupling  has 
been  developed,  but  still  has  some  features  that  are  not  well  understood.  The 
theory  predicts  “resonance  frequencies”  at  which  coupling  between  shear  and 
compressional  waves  may  be  greatly  enhanced.  The  resonance  frequencies  can 
be  fairly  high  (>1  kHz)  for  gradients  estimated  for  volume  inhomogeneities  in 
shallow  water  sediments.  The  effects  of  the  resonance  have  not  been  studied 
nor  has  a  stochastic  theory  of  gradient  driven  coupling  been  developed.  This 
report  concludes  with  some  recommendations  for  future  research  to  extend  these 
approaches  and  to  apply  them  to  understanding  scattering  from  the  subbottom. 


I.  INTRODUCTION 


Researchers  face  formidable  challenges  in  determining  the  major  acoustical 
processes  affecting  acoustic  propagation  in  shallow  water  and  in  developing  the 
databases  and  stochastic  models  needed  to  accurately  predict  the  acoustic  field. 
The  challenges  arise  from  the  complexity  of  the  acoustic  environment  (water 
column  and  subbottom).  The  environment  is  characterized  by  variability  on 
scales  that  affect  propagation  from  Very  Low  Frequencies  (5  Hz)  to  ultrasonic 
frequencies  (1  MHz)  frequencies.  Internal  waves,  fine  structure  of  the  water 
column,  subbottom  inhomogeneities  and  layers,  bottom  bathymetry  and 
roughness  have  all  been  suggested  as  significant  environmental  parameters. 

The  challenge  is  to  understand  the  physics  of  these  processes  well  enough  to 
develop  the  simplified  models  and  databases  needed  for  adequately  predicting 
the  performance  of  Navy  acoustic  systems. 

An  essential  component  of  the  solution  to  this  problem  is  an  understanding  of 
scattering  from  the  seafloor  and  subbottom.  Scattering  produces  the 
reverberation,  frequency  spread,  time  spread,  and  angle  spread  that  affect  naval 
acoustic  systems.  Yamamoto^  has  developed  a  stochastic  theory  of  scattering 
from  measured  velocity  and  density  inhomogeneities  that  produced  estimates  of 
backscatter  strength  that  are  consistent  with  measured  backscatter  in  the  1 00  Hz 
to  10  kHz  frequency  range.  Other  work  by  Anderson^  demonstrated  that  bottom 
backscatter  at  6.5  kHz  could  be  explained  by  scattering  from  density  and  velocity 
fluctuations  measured  in  cores.  The  fluctuations  in  density  and  velocity  in  the 
subbottom  in  shallow  water  are  important  sources  of  acoustic  scattering  that 
need  to  be  included  in  Navy  models  and  databases. 

To  develop  adequate  models  and  databases  for  predicting  these  phenomena,  an 
accurate  representation  of  seafloor  structure  (bathymetry,  bottom  roughness,  and 
subbottom  structure)  is  needed  on  the  scale  of  an  acoustic  wavelength.  Recent 
work  has  made  progress  toward  providing  structure  on  this  scale.  Yamamoto^ 
has  developed  the  technology  of  borehole  tomography  and  produced  the  first 
direct  measurements  of  velocity  inhomogeneities  (as  opposed  to  distinct  layers) 
in  the  upper  50  m  of  the  subbottom  on  scales  ranging  from  less  than  1  m  to 
about  10  m.  Goff  and  Jordan^  have  developed  a  stochastic  bathymetry  model 
that  has  the  potential  for  extrapolating  multibeam  bathymetry  data  to  millimeter 
scales.  Other  work  by  Orsi®  uses  x-ray  tomography  of  cores  to  provide  millimeter 
to  centimeter  scale  measurements  of  density  fluctuations. 

This  report  describes  research  carried  out  by  Science  Applications  International 
Corporation  from  1995  to  1998  to  develop  the  technology  for  understanding 
scattering  from  the  seafloor  and  subbottom  in  shallow  water  at  frequencies  from 
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100  Hz  to  10  kHz.  The  overall  objective  of  this  research  was  to  develop 
stochastic  acoustic  and  geological  models  for  accurate  numerical  simulation  of 
the  effects  of  scattering  from  the  seafloor.  There  were  three  technical  goals. 

First,  to  develop  a  quantitative  understanding  of  the  role  of  subbottom 
inhomogeneities  in  producing  reverberation  and  FAT  (frequency,  angle,  and 
time)  spreads  in  shallow  water.  Second,  to  develop  stochastic  models  describing 
subbottom  parameters  producing  scattering.  Third,  to  develop  acoustic  models 
for  predicting  FAT  spreads  and  reverberation  from  stochastic  geological  models. 
Due  to  a  funding  shortfall,  the  scope  of  our  work  was  reduced  and  focused  on  the 
first  goal. 

We  concentrated  on  two  theoretical  approaches  that  looked  promising  for 
application  to  scattering  from  subbottom  inhomogeneities.  The  first  was  the  T- 
matrix  approach  to  scattering.  We  examined  this  technique  and  found  that  it  has 
all  the  attributes  necessary  for  application  to  scattering  from  objects  buried  in  the 
subbottom.  The  second  was  the  coupled  potential  approach  for  predicting 
gradient  driven  coupling  between  shear  and  compressional  waves.  This 
approach  appears  suitable  for  application  of  the  Born  approximation  for 
predicting  scattering  from  subbottom  velocity  and  density  fluctuations.  Both  of 
these  approaches  are  reviewed  below. 

As  a  result  of  our  work,  we  recommend  that  additional  research  be  carried  out  to 
develop  the  numerical  implementation  of  both  of  these  approaches  for  the 
prediction  of  time  and  angle  spreading  of  acoustic  signals. 

II.  PRELIMINARY  RESEARCH 

We  began  our  research  with  some  two  preconceived  notions.  The  first  of  these 
was  that  we  were  looking  for  approaches  with  the  fewest  theoretical  limitations. 
Our  goal  was  to  develop  an  approach  that  could  be  used  for  sensitivity  studies  to 
identify  the  major  mechanisms  and  parameters  that  would  be  needed  for 
modeling  the  performance  of  naval  systems.  To  carry  out  these  studies  we  felt 
that  accuracy  was  paramount  in  the  research  model  so  that  as  many  physical 
mechanisms  as  possible  were  included.  The  sensitivity  studies,  not  our  model, 
should  identify  the  important  mechanisms  and  parameters.  The  second  notion 
was  that  the  large  gradients  associated  with  volume  inhomogeneities  might  make 
gradient  driven  coupling  of  shear  and  compressional  waves  a  significant  process 
in  shallow  water.  This  idea  was  based  on  the  recent  reports  that  velocity 
fluctuations  could  have  gradients  an  order  of  magnitude  larger  than  those 
previously  considered. 
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III.  T-MATRIX  APPROACH  TO  SCATTERING 


Our  research  began  with  a  survey  of  recent  publications  in  the  Journal  of  the 
Acoustical  Society  of  America  to  identify  potential  approaches  to  scattering.  By 
far,  we  found  that  the  majority  of  publications  between  1976  and  1994  dealing 
with  scattering  were  based  on  the  T-matrix  approach  or  some  variant  of  it.  We 
found  publications  dealing  with  scattering  of  acoustic  waves,  scattering  of  elastic 
waves,  multiple  scattering,  scattering  by  various  shaped  objects,  scattering  from 
objects  with  internal  resonances,  and  scattering  from  objects  in  a  waveguide. 
While  the  T-matrix  approach  has  been  applied  to  scattering  from  objects  buried 
in  the  seafloor,  it  has  not  yet  been  applied  to  scattering  from  sedimentary 
layering  in  the  subbottom.  It  appeared  from  this  survey  that  the  T-matrix 
approach  should  be  examined  for  its  suitability  for  predicting  scattering  from 
subbottom  inhomogeneities. 

Waterman®  originally  developed  the  T-matrix  approach  for  acoustic  scattering 
problems.  The  approach  uses  the  Helmholtz  integral  equations  that  relate  the 
scattered  field  to  the  incident  field  and  an  integral  of  the  total  field  on  the  surface 
of  the  scattering  object.  The  incident  and  scattered  fields  are  expanded  in  terms 
of  a  set  of  basis  functions,  such  as  spherical  harmonics.  The  essence  of  the  T- 
matrix  approach  is  to  use  two  surface  integral  equations,  one  for  an  observation 
point  inside  the  scattering  object  and  the  other  for  a  point  outside  the  scattering 
object,  to  obtain  a  relationship  between  the  expansion  coefficients  of  incident  and 
scattered  fields.  The  matrix  relating  the  coefficients  is  the  transition  matrix,  the 
“T-matrix”.  The  T-matrix  depends  only  on  the  shape  of  the  scattering  object  and 
the  boundary  conditions.  The  only  mathematical  limitations  imposed  by  the 
theory  are  that  the  shape  of  the  scatterer  allow  an  inscribed  and  circumscribed 
sphere  for  defining  the  spherical  wave  functions  used  as  a  basis  set  for  the 
incident  and  scattered  field.  Otherwise,  the  scattering  surface  can  have  almost 
any  shape. 

The  T-matrix  approach  proved  to  be  very  general  and  has  been  applied  to 
scattering  of  elastic  and  electromagnetic  waves.  In  these  cases,  the 
mathematics  is  more  complicated  than  that  for  acoustic  scattering  because  one 
must  deal  with  vector  fields  and  tensor  Green’s  functions.  Ref.  7  is  a 
comprehensive  listing  of  publications,  as  of  1988,  extending  and  applying  the  T- 
matrix  approach  and  Ref.  8  is  collection  of  early  work  on  the  T-matrix.  A  unified 
view  of  T-matrix  theory  for  acoustic,  elastic,  and  electromagnetic  problems  is 
given  in  Ref.  9  while  a  recent  review  article^®  concentrates  on  application  to 
scattering  from  structured  objects  with  resonances. 
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For  the  research  reported  here,  we  examined  papers  on  the  application  of  T- 
matrix  theory  to  elastic  scattering®’^ scattering  from  objects  with 
attenuation^^’^®,  multiple  scattering^®'^®,  scattering  from  structured  objects^ 
scattering  from  non-spherical  objects^^'®^,  scattering  from  objects  in  a  layered 
medium®®'®^,  and  scattering  in  a  waveguide®®'®^.  Our  major  focus  was  to 
understand  the  theory  well  enough  to  determine  the  feasibility  of  applying  it  to 
scattering  from  subbottom  inhomogeneities.  We  did  not  identify  any  fundamental 
limitations  that  would  exclude  it  from  this  application. 

A.  Acoustic  scattering 

We  will  follow  the  approach  of  Waterman®  in  briefly  outlining  the  derivation  of  the 
fundamental  T-matrix  equations.  Waterman  examines  the  total  acoustic  field, 
being  made  up  of  an  incident  field,  and  a  scattered  field,  'P®.  His  starting  point 
is  the  Green’s  theorem  for  scalar  fields®  relating  the  fields  inside  and  outside  a 
bounding  surface  to  the  sources  of  the  field  and  the  field  on  the  surface.  This  is 
essentially  Huygen’s  principle.  The  geometry  is  shown  in  Fig.  1.  There  are  no 


sources  within  the  scatterer.  In  this  case  the  bounding  surface  is  that  of  the 
scattering  object  (o)  and  a  sphere  at  infinity.  The  volume  we  consider  is  all 
space  excluding  the  interior  of  the  scattering  object,  i.  e.,  the  interior  of  the 
scatterer  is  the  exterior  of  our  volume.  On  the  surface  (a)  of  the  scatterer,  the 
outward  normal  of  the  surface  bounding  our  volume  is  ih  which  is  in  the  opposite 
direction  to  the  outward  normal  of  the  surface  bounding  the  scatterer,  n .  The 
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source  is  has  a  single  frequency  and  has  an  exp(-/(jot)  time  dependence.  There 
are  two  integral  equations,  one  for  an  observation  point  in  our  volume  and  one 
for  an  observation  point,  outside  our  volume.  For  a  point  within  the  volume. 
Green’s  theorem  is 

W  =  '¥‘(r)+-^jdan- ['F^(r')Vg( *,|r - r'|)  - *,|r - r1) V^'K(r’)]  (1) 


where:  r  is  an  observation  point  inside  the  volume,  r'  is  a  point  on  the  surface  a, 
kp  =  (ojcp  is  the  wavenumber  of  the  acoustic  field,  is  the  sound  speed  of  the 
medium,  and  V^.'F  are  the  total  field  and  its  gradient  on  the  surface  of  the 
scatterer  (on  the  outside  of  the  scatterer),  and  g  is  the  free  space  scalar  Green’s 
function.  The  free  space  Green’s  function  satisfies  V^g  +  A:^g  =  4;r<5(r-r')  and  is 
given  by 


ikpR 

i(k,R)=ik,hS^(/,,R)=^ 

where:  /?  =  |r  -  r'|  and  /j^'^  is  the  Hankel  function  of  the  first  kind  of  order  zero. 
The  Green’s  function  represents  an  outgoing  spherical  wave  from  a  source  at  r' . 

The  second  integral  equation  is 

0=  -  jr(i,|r-r'|)V^>P(r')]  (2) 


where  the  observation  point  r  is  now  outside  the  volume  (inside  the  scatterer). 
This  equation  states  that  the  incident  field  and  the  field  generated  by  the  virtual 
sources  on  the  surface  cancel  inside  the  volume.  It  does  not  say  that  the  total 
field  inside  the  volume  is  zero.  The  field  in  the  interior  is  given  by  the  solution  to 
the  wave  equation  for  the  interior  with  boundary  conditions  and  total  field 
specified  on  the  surface  of  the  scatterer. 

The  next  step  is  to  expand  the  Green’s  function,  the  incident  field,  and  the 
scattered  field  in  terms  of  basis  sets  whose  coefficients  can  be  related  through 
the  boundary  conditions  on  the  surface  of  the  scatterer.  Waterman  chose  the 
partial  wave  solutions  of  the  scalar  Helmholtz  equation,  V^'F  + =  0 ,  in 
spherical  coordinates.  (He  also  considered  cylindrical  coordinates.)  In  this  case, 
the  basis  set  is  given  by 
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where:  /j^'^  are  spherical  Hankel  functions  of  the  first  kind  of  order  n,  are 
spherical  harmonics,  p  is  an  index  with  two  values  (even  or  odd)  indicating  the 
parity  of  the  spherical  harmonic,  and  „  is  a  constant  given  by 

f  {in  + 1)  for  w  =  ol 

|2(2w  +  lX»  -  w)!/(«  +  m)\  for  m  Oj 

The  spherical  harmonics  are  given  by 

I-.'. (eA)= 

[sin  {m<f)P”  (cos  0)  for  p  =  odd  J 


where  PJ”  are  the  associated  Legendre  polynomials.  We  choose  a  coordinate 
system  with  an  origin  inside  the  scattering  object.  Since  there  are  no  sources 
within  the  scattering  object,  we  expand  the  incident  field  in  terms  of  regular  wave 
functions  i.,  e.  wave  functions  that  do  not  have  singularities  at  the 

origin.  For  spherical  coordinates  the  regular  wave  functions  are 

Re  =  y'mjn  ^  pr)^ l„  {O, 


The  incident  field  is  then 


'J"  =I]««ReV'„ 


in  which  we  have  collapsed  the  multiple  indices  (p,  m,  n)  to  a  single  index,  n. 
Expansions  for  plane  waves  and  point  sources  in  homogeneous  media  are 
available  for  spherical  and  cylindrical  coordinate  systems. 
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The  Green’s  function  for  free  space  may  be  expanded  in  terms  of  the 
eigenfunctions  as^ 


S^P  l«*  -  *'l) = X  (^p^>  ) 

where:  and  refer  to  the  lesser  and  greater  of  r  and  r'.  If  we  consider  the 

region  outside  a  sphere  enclosing  the  scattering  object,  express  the  scattered 
field  as  the  expansion 


use  the  eigenfunction  expansion  for  the  Green’s  function,  and  note  that  r’>r,  then 
Eq.  (1)  becomes 

/.  =  ^jrfcTn  •[4'^(r')VRe»-,(*,r')-  (3) 

Similarly,  for  a  region  inside  a  sphere  contained  within  the  scattering  object,  Eq. 
(2)  becomes 


To  proceed  further,  we  need  to  consider  the  boundary  conditions  on  the 
scattering  object.  The  boundary  conditions  relate  the  total  field  to  its  gradient  on 
the  surface  of  the  object.  We  will  consider  two  boundary  conditions:  the 
Neumann  boundary  condition  (n •  v  T'  =  0 ),  and  a  penetrable  boundary  condition 
having  continuous  pressure  and  normal  component  of  velocity.  Expanding  the 
total  field  in  a  complete  set  of  functions  allows  us  to  solve  Eqs  3  and  4  can  then 
for  a„  and  /„.  The  matrix  relating  them  is  the  T-matrix. 

For  the  Neumann  boundary  condition,  we  choose  to  expand  the  field  on  the 
surface  of  the  object  in  terms  of  the  regular  eigenfunctions,  which  we  assume 
forms  a  complete  set  on  the  surface  of  the  object,  as 
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Substituting  into  Eqs.  3  and  4  yields  the  following  two  equations  in  matrix  and 
vector  notation. 


a  =  iQ^a  (5) 
f  =  /Re(Q^)a  (6) 

The  matrix  Q  is  given  by 

Q”-"  "  ^ J 

ReQ^„  =  ^  Jfifo-n  •v(Re  y/^  (t^r '))Re  if/„  [k/) 

Eliminating  a  from  Eqs.  5  and  6  yields  the  T  matrix 

f  =  Ta 


where 

T  =  -Re(Q^XQ")'' 


The  next  example  is  the  usual  acoustic  boundary  conditions  on  the  surface  of  the 
scattering  object— continuity  of  pressure  and  the  normal  component  of  the 
velocity.  Primes  will  denote  quantities  in  the  interior  of  the  object.  The  boundary 
conditions  on  the  surface  are 

n.V,^  =  n.V_'R 

Assuming  the  interior  to  be  homogeneous,  the  wavefunctions  Jle'P„(yt'r)}form  a 
complete  set  and  the  field  in  the  interior  can  be  expanded  as 
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Using  the  boundary  conditions  to  express  in  terms  of  'F_,  the  eigenfunction 
expansion  of  and  eliminating  the  coefficient  vector  p  results  in  the  following 
expression  for  the  Q  matrix  that  defines  the  T-matrix. 

Attenuation  can  be  included  by  making  A:'  complex. 

The  final  step  is  the  evaluation  of  the  infinite  size  T-matrix.  This  can  be  done  by 
solving  for  the  T  matrix  of  rank  n  and  then  increasing  n  until  the  matrix  elements 
converge.  Reference  10  discusses  more  numerically  efficient  approaches. 

The  advantages  of  the  T-matrix  over  other  approaches  to  scattering  can  be  seen 
from  the  examples  above.  The  equations  forQ,  and  hence  the  T-matrix  depend 
only  on  the  size  and  shape  of  the  scatterer  and  not  on  the  incident  field.  Once 
evaluated,  the  T-matrix  can  be  used  to  find  the  scattered  field  produced  by  any 
incoming  field.  Another  advantage  is  that  the  methodology  is  not  restricted  to 
objects  of  a  certain  shape.  The  integral  involved  is  a  surface  integral.  Once  an 
integration  algorithm  is  developed,  it  can  be  applied  to  any  surface.  Unlike  the 
eigenfunction  expansion  method,  the  T-matrix  is  not  restricted  to  objects  of  a 
particular  shape  having  the  symmetry  of  the  coordinate  system. 

B.  Elastic  wave  scattering 

The  T-matrix  for  elastic  wave  scattering  is  much  more  complicated  than  that  for 
acoustic  scattering.  The  basic  reason  for  this  is  that  the  basis  functions  become 
a  set  of  vectors,  rather  than  a  single  scalar,  and  the  Green’s  function  becomes  a 
tensor  quantity.  For  our  discussion  below  we  follow  the  presentation  in  Ref.  14. 

We  consider  an  infinite,  homogeneous,  isotropic  medium  of  density  p  and  Lame 
parameters  X  and  p.  The  embedded  scatterer  has  different  properties  pi,  Xi,  and 
Pi.  The  only  condition  on  the  surface  of  the  scatterer  is  that  it  be  smooth,  having 
a  normal  vector  that  changes  continuously  over  the  surface.  The  geometry  is  still 
that  of  Fig.  1 .  The  incident  wave  is  time  harmonic  with  frequency  ©  and  time 
dependence  The  vector  displacement  field  u(r)  is  made  up  of  the  incident 
displacement  u'(r)  and  the  scattered  displacement  u^(r)  where  r  is  the  position 
vector  of  an  observation  point.  The  equation  of  motion  is 

V-’r+p<y^u  =  0  (7) 
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where  T  is  the  stress  tensor,  which  is  defined  in  terms  of  the  displacement  as 

T  =  XiV  •  u  +  //(Vu  +  uV) 


i  is  the  identity  tensor.  The  traction  t  is  given  by  t  =  n.  •  f . 

The  first  step  in  the  derivation  of  the  T-matrix  for  elastic  wave  scattering  is  to 
obtain  integral  equations  similar  to  Eqs.  1  and  2.  These  integral  equations  relate 
the  incident  and  scattered  field  to  the  integral  over  the  surface  of  the  scatterer. 
We  begin  by  defining  the  Green’s  tensor  G  and  a  Green’s  stress  tensor  S  as 
being  related  to  a  point  source  at  r'.  The  Green’s  tensor  is  the  solution  to 

VS  +  p(y'G  =  -i(5(r-r')  (8) 


where  S  is  the  stress  tensor  (rank  3)  defined  in  terms  of  the  Green’s  tensor 


S  =  AiV-G  +  //(vG  +  Gv) 


In  indicial  notation,  the  last  term  in  the  definition  of  S  is  defined  as  +djGi^ 
where  5,  =d/dx.  and  /,  J,  and  k  refer  to  the  components  of  the  position  vector. 
The  solution  of  these  two  equations  is  the  Green’s  tensor  given  by®'^® 

-  *'1)-  [g(k^\r-r])- g(ArJr  -  r'|)] 


Where:  g  is  the  scalar  Green’s  function  defined  above,  =  p(o^/{A  +  2ju)  is  the 
square  of  the  compressional  wavenumber,  and  jp  is  the  square  of  the 

shear  wavenumber. 

The  integral  equation  relating  the  displacement  field  to  the  field  on  the  surface  of 
the  scatterer  is  obtained  from  Eqs.  7  and  8.  We  apply  the  divergence  theorem  to 
Eq.  7  postmultiplied  by  G  and  Eq.  8  premultiplied  by  u  to  obtain 

u(r)=u'(r)+ Jcfcrlu-n-S-n-T-G}  for  r  outside  <7  (9) 

0  =  u'(r)+ ffi?cr{u'n  S-n-T-G)  for r inside  cr  (10) 
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We  note  that  the  quantities  inside  the  integrals  depend  on  both  r  and  r'  where 
r'  is  a  point  on  the  surface  and  that  da  is  the  area  element  at  r'. 

The  next  step  in  obtaining  the  T-matrix  is  to  define  the  basis  functions  we  will  use 
to  expand  the  field  quantities.  In  three  dimensions  the  displacement  vector  can 
be  written  as  the  sum  of  three  vectors  derived  from  scalar  potentials  as 

u(r)  =  VP  +  ky  X  (r2)+  V  x  V  x  (r5) 


where:  compressional  waves  are  described  by  the  potential  P  and  shear  waves 
by  the  potentials  Q  and  5 .  In  analogy  to  this  equation  there  are  three  sets  of 
basis  functions 


^  A.  (cosP)cos»i97] 

.(^p  l^s  Y'  cMK  (^pi'h"  (cos  0)cosm<p] 


p  =  even 
p  =  odd 


X  [rh„{kj)P„"'  (cos0)cos/n^] 
k^7j„V  X  [r/2„(A:^r)P”(cos^)cosm^] 


p  =  even 
p  =  odd 


Xpnm(^)  =  -r^^¥pnn,{^) 


Arrows  over  quantities  in  the  above  equations  indicate  that  they  are  vectors.  The 
constants  and  are  give  by 


(2«  +  lX«-/w)! 
4;r(«  +  /w)! 


{in  +  lX«  -  ni)\ 

"  Amin  +  iXw  +  m)\ 


where  =  1  and  £■„  =  2  (n  >  0)  is  the  Neumann  factor.  Basis  functions  that  are 
regular  at  the  origin  are  obtained  by  replacing  h„  with  . 
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The  quantities  appearing  in  Eqs.  9  and  10  can  now  be  expanded  in  terms  of  the 
basis  functions. 

pnm  Re??,„„(r)+C^„„  Ref  ^„„(r)] 

pnm 

ocpnn,  Re^p„„(r)+>9,„„  Re<?p„„(r)+^  Ref  ;,„„(r) 

pnm 

*■ ')  =  7^-V  X  fc”™  ('*<  ^P"'”  (*■> ^p»’”  (*■<  (*■>  )+X  pnm  (»•<  )Re  f  (r>  )] 

\pf^)  ^ 


The  equation  for  G(r,r')  follows  from  the  expansion  of  the  scalar  Green’s 
functions  in  terms  of  the  basis  sets  in  the  expression  for  G,^.  given  above. 

The  integral  equation  for  points  inside  a  sphere  within  the  scattering  volume  (Eq. 
9)  becomes,  noting  that  r  <r’ , 

pnm  Re^p„„(r)+^p„„  Rev^;,„„(r)+Cp„„  Ref  ^„„(r)]=  (1 1) 

pnm 

"  2  J  n  •  r')-  n  •  t(r')-  G  ^„„(r,  r')} 

pnm 


where  the  subscripts  on  the  Green’s  tensor  and  Green’s  stress  tensor  denote 
terms  in  their  expansions  for  r  <  r'  given  below. 

G('‘>  •*')  =  7-^  ^  (•■  ')Re  hnm  («*)  +  ¥  pnm  ')Re  ¥  pnm  fr)  +  X  pnm  (rORe  X  pnm  (»•)] 

\P^ )  pnm 

+  ^  ^P’'”'  ('■')^']Re  ¥  pnm  (r)  +  /^[V'f  p„„  (r ')  +  Xp„m  ']Re  Xp„m  (r) 
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Applying  the  same  procedure  to  Eq.  10,  and  recalling  that  r  <r'  for  this  case, 
yields  the  following  equation  for  the  expansion  coefficients  of  the  scattered  field. 

X  («*)  +  P pnJ' pnm  (»*)  +  Y  pnmX pnm  fr)]  =  (12) 

pnm 

-  X  ^p'"” 

pnm 


where  G^„„(r',r)  is  obtained  by  interchanging  r  and  r'  in  the  above  expression 
for  G^„„(r,r')  and 

(po)) 

+ T^r  (■■') + (■■) 

+  ^[V' Re  f  ^(r')+ Re  ^  (r')V'l?,..  (r)} 


Equations  11  and  12  do  not  reduce  to  the  simple  relationship  between  an 
expansion  coefficient  and  a  surface  integral,  as  was  the  case  for  acoustic 
scattering  (Eqs.  5  and  6).  To  obtain  a  similar  result  for  elastic  scattering,  we 
multiply  Eq.  11  successively  by  Re^^„„(r),  Re5?p„„(r),  and  Ref^„„(r),  integrate 
over  all  angles,  use  the  orthogonality  relations  of  spherical  harmonics,  and 
perform  a  little  algebra  to  obtain 

4-  =  (*“')+  AVVp„„(r')+  A^;,„„(r')V'] 

-n-t(r')-^,„„(r')} 

J* dG^{r')- n  •  [/iV V^„„(r')+  n  •  t(r')- 

MXpnS^’y]- n  •  T(r')- 


Equation  12  yields  similar  expressions,  with  the  basis  functions  replaced  by 
regular  basis  functions,  for  the  coefficients  of  the  scattered  field. 
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As  was  the  case  for  acoustic  scattering,  we  need  to  employ  the  boundary 
conditions  on  the  scatterer  to  proceed  further.  As  an  example,  we  consider  a 
cavity  with  stress  vanishing  on  the  surface  of  the  scatterer.  The  boundary 
condition  is  n-f(r')=0  for  r'  on  o.  We  expand  the  surface  displacement  in 
terms  of  the  regular  basis  set  as 


in  which  we  have  used  different  subscripts  to  identify  the  basis  functions  and 
coefficients.  After  substituting  this  equation  into  the  expressions  for  the 
coefficients  of  the  incident  and  scattered  fields  we  obtain  two  matrix  equations. 
One  relating  the  incident  wave  coefficients  to  the  coefficients  of  the  surface  field 
and  the  other  relating  the  scattered  wave  coefficients  to  the  surface  field 
coefficients. 

fe't  (e'‘t 

fc'KL  (e”)r.. 


A,..  =-'Re(e"')l  Refect  K 

r„.\  [Refc”);!  Refect 

Each  vector  component  in  these  equations  is  an  infinite  component  vector  itself. 
Each  matrix  element  is  an  infinite  dimension  matrix.  The  for  fh®  cavity 
problem  are  given  by  these  integrals  over  the  surface  of  the  scattering  object 

fe"  t  (■•')■ « •  fllV'  • 

fe”  t  =  ‘'<^Re4(>-’)'  A  •  [V'f  ,™,(r')+ 


The  matrix  elements  are  obtained  by  replacing  Re^^,^(r')  with  Re^^.(r') 
and  (|2*")  are  obtained  by  replacing  Re^,^(r')  with  Re^^(r').  Re((2pL)  is 
obtained  by  replacing  the  basis  vectors  by  the  corresponding  regular  basis 
vectors. 

The  T-matrix  is  obtained  by  solving  the  two  matrix  equations  above  for  the 
coefficients  of  the  scattered  wave  in  terms  of  the  coefficients  of  the  incident  ware. 


a 


pnm 


pnm 


pnm 


/pnm 


(T-’t 


The  elements  of  the  T-matrix  are  defined  in  terms  ofQby  the  same  equation  as 
in  the  scalar  case 


T  =  -Re(Q^XQ")'' 


As  in  the  acoustic  scattering  case,  the  T-matrix  depends  only  on  the  geometry 
and  boundary  conditions  of  the  scatterer,  not  on  the  properties  of  the  incident 
field.  More  complicated  boundary  conditions,  such  as  those  describing  an  elastic 
scatterer,  are  also  discussed  in  Ref.  14. 

C.  Multiple  Scattering 

The  T-matrix  approach  lends  itself  to  describing  multiple  scattering—the 
interaction  of  the  scattered  field  with  another  obstacle.  The  key  to  treating 
multiple  scattering  is  the  ability  to  translate  the  origin  of  coordinates  of  the 
spherical  basis  functions.  This  allows  one  to  calculate  the  scattering  from  each 
individual  object  as  outlined  above  and  then  to  translate  the  results  to  a  common 
coordinate  system.  For  multiple  scattering,  the  scattered  field  of  one  object 
becomes  the  incident  field  at  another.  In  this  section,  we  follow  Ref.  19  in 
discussing  the  application  of  the  T-matrix  to  multiple  scattering. 

1 .  Two  Scatterers 

Figure  2  shows  the  geometry  for  two  scatterers  Si  and  S2.  There  is  an  origin  of 
coordinates  0  outside  the  scatterers  and  points  Oi  and  O2,  located  by  the 
vectors  a,  and  a^that  serve  as  local  origins  of  coordinate  systems  within  the  two 
scatterers.  The  vectors  r,'  and  are  points  on  the  surfaces  and  02  of  the 
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scatterers  relative  to  the  O  coordinate  system  while  r,"  and  r"  are  points  on  the 
surfaces  relative  to  the  coordinate  systems  within  each  scatterer.  The  vectors  r 
are  observation  points  relative  to  the  0  system  while  r,  and  are  observation 
points  relative  to  the  Oi  and  O2  coordinate  systems.  Oi  and  O2  are  chosen  so 
that  r,"  and  r"  are  continuous  functions  of  spherical  coordinates  in  their 
respective  coordinate  systems.  The  vectors  lij  and  n2  are  outward  normal 
vectors  at  the  surface  of  each  scatterer. 


In  analogy  with  the  single  scatterer,  we  have  two  integral  equations  for  the 
acoustic  field  one  for  r  outside  Si  and  S2 

'P'(r)  =  ^  pOT-[T^(r')Vj(i>-r1)-*(i,|r-r1)V^T(rO]  (13) 

5|+5’2 

and  the  other  for  r  inside  Si  or  S2 

(r)  =  ^dan* (r ')Vg( |r - r'|)  - |r - r '|) V^^(r')]  (14) 

^1+^2 


In  Eqs.  1 3  and  1 4  we  assumed  that  there  are  no  sources  inside  a  sphere 
centered  at  0  and  enclosing  both  scatterers.  As  for  the  acoustic  scatterer  above, 
we  expand  the  incident  field  inside  this  sphere  as 


'P'(r)-^a„Rev^„(r)  (15) 


16 


To  proceed  we  will  want  to  expand  the  surface  fields  in  Eqs.  13  and  14  in  terms 
of  regular  eigenfunctions  centered  at  Oi  and  O2.  We  will  also  need  to  express 
the  expansion  for  in  terms  of  coordinates  centered  at  Oi  and  O2.  As  an 
example,  consider  the  case  for  r  inside  Si.  We  use  the  fact  that  r-r/  =  rj  -r* 
to  expand  the  Green’s  function  within  a  sphere  within  Si  centered  at  Oi.  We  also 
restrict  r  to  those  points  inside  Si  for  which  </'J^  where  is  the  smallest 
value  of  r"  on  the  surface  ai.  We  can  now  identify  and  in  the  expansion  of 
the  Green’s  function  to  obtain 


We  also  need  to  expand  the  Green’s  function  in  the  integral  over  02.  This  is 
accomplished  by  noting  that  the  argument  of  this  Green’s  function  (r-rj)  can  be 
written  as  r, -[r;-(a, -a,)].  For  /•,  <r,l,  we  have  <|r;-(a, -aj| .  The 
Green’s  function  is 


s(kp  |r  -  r'  1)  =  ik^  2]  Vp  («•;  -  (a,  -  a,  ))]Re  if/„  {k^r, )  (17) 


To  carry  out  the  integral  over  oi  we  need  the  translation  properties  of  y/^  and 
Re^„  to  convert  the  arguments  to  the  variables  of  integration.  For  ,  noting 
that  r  =  r,  +a, ,  the  translation  is  given  by  a  summation  over  regular  functions, 

Rey/^Cr,  +ai)=^7?„„,(a,)Rey/„,(r,)  (18) 


The  translation  properties  of  y/^  is  also  given  in  terms  of  the  regular 
eigenfunctions  as 


+a2)Rey^„,(r3")  (19) 

rt* 

This  expansion  is  valid  only  for  [a,  -a2|  <^2" .  For  the  expansion  inside  Si  to  be 
valid  we  also  require  |aj  -a2|<A',''.  These  two  conditions  introduce  a  constraint 
on  the  configuration  of  scatterers  that  can  be  satisfied  if  the  spheres  defined  by 
r,"  and  do  not  overlap.  The  properties  of  the  translation  matrices  and 
r„„,  are  given  in  Ref.  19. 
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Equations  18  and  19  can  now  be  used  in  Eqs.  13  and  14  to  relate  the  incident 
and  scattered  fields  to  the  fields  on  the  surface  of  the  scatterers.  The  expansion 
coefficients  of  the  surface  fields  are  a‘  where  /  =  1,2  denotes  the  scattering 
object.  The  Q  matrices  for  each  scattering  object,  as  defined  above  for  acoustic 
scattering,  will  also  appear  in  the  results.  A  superscript  will  identify  the  matrix  for 
each  scattering  object  as  Q'  .  The  supersaipt  T  denotes  the  transpose  of  a 
matrix.  Equations  13  and  14  then  yield  the  following  equations  in  matrix  notation. 
For  r  =  r,  +a,  inside  Si  with  <r,„^ , 


R^(a2)a  =  iQ'a'  +/T(-a,  +a2)ReQ^a^  (20) 
and  for  r  =  +  a^  inside  Si  with  . 

R^(a2)a  =  /r(-a2  +a,)ReQ'a'  +iQ^a^  (21) 

For  the  region  outside  the  scatterer  we  obtain  an  equation  for  the  scattered  field 
expansion  coefficients. 

f  = -/R(ai)ReQ'a'  -iR(a2)ReQ^a^  (22) 


Equations  19-21  can  be  used  to  eliminate  the  expansion  coefficients  of  the  field 
on  the  surface  of  the  scatterers  and  obtain  the  T-matrix  for  the  two  scatterers 
t(i,2).  The  T-matrices  of  the  individual  scatterers,  T(/)=-ReQ*(Q'  ) ',  also 
appear  in  expression  for  the  total  T-matrix.  The  T-matrix  for  two  scatters  is 

t(1,2)=  R(a,){r(lXl-r(-.,  +a,)r(2X'(-a,  +a,)r(l)r  (23) 

x[l  +  r(-a,  +ii^)T{2)R(a,  -a2)]R(-a,)} 

+  R(a>){r(2Xl  -r(-a,  +a,)r(l)r(-a,  +a,)r(2)]-' 
x[l  +  r(-a,  +a,)r(l)R(a,  -a,)]R(-a,)) 

Examination  of  Eq.  23  shows  that  it  treats  the  two  scatterers  in  a  symmetric  way. 
The  introduction  of  more  than  one  scatterer  does  not  require  more  complicated 
surface  integrals  than  the  single  scatterer  case.  We  still  need  only  calculate 
theQ' matrices  of  the  individual  scatterers.  Equation  23  combines  them  into  the 
T-matrix  for  2  scatterers. 
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It  is  informative  to  examine  some  of  the  individual  terms  arising  from  Eq.  23.  For 
T(i)->0  we  recover  the  single  scatter  result  T(l,2)->R(a,)T(l)R(-a])  which  is 
the  single  scatter  result  translated  from  an  origin  inside  the  scattering  object  to 
one  outside  the  object.  Formally  expanding  the  inverses  in  Eq.  23  yields  a  series 
that  has  a  natural  interpretation  in  terms  of  multiple  scattering.  The  first  few  of 
these  terms  are  interpreted  in  Table  1 . 


Expansion  term 

interpretation 

T(1) 

Single  scatter  from  1 

_  T(2) 

Single  scatter  from  2 

T(2)rT(l) 

Scatter  from  1 ,  propagation  from  1  to  2, 
scatter  from  2 

T(1)rT(2) 

Scatter  from  2,  propagation  from  2  to  1 , 
scatter  from  1 

T(i)rT(2)rT(i) 

Scatter  from  1 ,  propagation  from  1  to  2, 
scatter  from  2,  propagation  from  2  to  1 , 
scatter  from  1 

T(2)rT(1)rT(2) 

Scatter  from  2,  propagation  from  2  to  1 , 
scatter  from  1 ,  propagation  from  1  to  2, 
scatter  from  2 

Table  1 .  Interpretation  of  terms  in  the  expansion  of  T(1 ,2) 


This  approach  to  multiple  scattering  is  very  powerful.  It  is  a  full  physics  approach 
that  includes  effects  of  energy  trapping  and  resonances  in  the  configuration  of 
the  scatterers. 

2.  N  Scatterers 

While  the  approach  for  2  scatterers  that  we  used  above  can  be  generalized  to 
treat  N  scatterers,  the  resulting  iterative  procedures  outlined  in  Refs.  19  and  22 
are  complex.  Reference  23  provides  a  different  approach  that  is  uses  the 
"effective  field”  concept  to  develop  a  concise  and  elegant  solution  for  the  T-matrix 
for  N  scatterers.  The  incident  field  at  each  scatterer  is  expanded  as 
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The  scattered  field  is  expanded  as 


'J"W=  l]'P;(r-<l;)=  V.' 

7=1  7=1  n 


We  define  the  effective  field  incident  on  each  scatterer  as  the  field  made  up  of 
the  incident  field  and  the  fields  scattered  from  all  other  obstacles.  Since  this  is 
the  field  exterior  to  the  scatterer  under  consideration,  it  can  be  expanded  in  terms 
of  regular  functions  centered  on  the  scatterer.  The  effective  field  incident  upon 
the  j"'  scatterer  is  then 

>Ff(r-d,)=r(r-d,)+^'P'(r-d,) 

j*J 

=  y'Re^/„(r-d^.)c;; 


The  expansion  coefficients  fj  of  the  field  scattered  by  the  f'  object  is  given  by 
the  T-matrix  for  that  scatterer  as 

<24) 


Expanding  the  expression  for  the  effective  field  and  using  the  translation  matrices 
defined  above  yields 


7*7 


In  matrix  notation  this  equation  can  be  written  as 


c  =  a  +  (25) 
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where 


C'n 

yO) 

0  ... 

I - 

O 

cl 

I! 

II 

0 

... 

0 

\ - 

\ _ 

1 

1 

O 

0  ... 

- 1 

H 

and 


0 

r(-(<i,-d,))  r(-(d,-d,))  ... 

r(-(d,  -d„))' 

r(d,-d,) 

0 

r(-(<i,-d,))  ... 

r(-(d,-d„)) 

r(d,-d,) 

r(d,-d,) 

0 

r(-(d,-d„)) 

1 

1  : 

a 

r(<i,-d„) 

r(d,-d„)  ... 

0 

Using  the  matrix  analog  of  Eq.  24,  c  =  7'  'f ,  in  Eq.  25  yields  the  multiple-scatter 
T-matrix  relating  the  incident  and  scattered  fields 

f  =  T^a 


where 


Expanding  the  inverse  in  produces  terms  corresponding  to  increasing  orders 
of  multiple  scattering 


Each  7'  matrix  represents  scattering  by  each  of  the  N  obstacles  while  each 
matrix  distributes  the  scattered  fields  among  the  N  objects.  The  poles  of 
correspond  to  scattering  resonances  of  the  N  body  system. 
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D.  Layered  Scatterer 

The  T-matrix  approach  has  been  extended  to  treat  scattering  from  a  structured 
object,  such  as  an  object  made  up  of  layers  of  different  materials,  an  object  with 
embedded  scattering  objects,  or  an  object  made  up  of  non-enclosing  parts.  The 
theoretical  approach  we  present  here  treats  scattering  from  a  layered  object  and 
follows  the  work  in  Ref.  24.  Figure  3  shows  the  scattering  geometry.  There  is  a 
single  object  with  two  layers.  A  subscript  / ,  /  =  1,2 ,  will  be  used  to  identify 
quantities  in  the  two  layers  of  the  scatterer.  A  subscript  0  will  denote  the 
properties  of  the  background  medium  outside  the  scatterers. 


Figure  3.  Geometry  for  a  layered  scatterer. 


We  begin  by  applying  Huygen’s  principle  to  the  region  outside  region  1  and  using 
the  acoustic  boundary  conditions  of  continuity  of  pressure  and  velocity  to  relate 
the  field  in  the  exterior  (Region  0)  to  the  interior  (Region  1).  The  incident  and 
scattered  fields  are  expressed  as 

(*«'■) 


n 
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In  the  region  outside  02  and  inside  oi  the  field  is  made  up  of  both  regular  (Re) 
and  outgoing  (Out)  solutions  and  is  written  as 


'F’(A:,r)  =  SI 

n 


The  integral  equations  yield  two  equations  relating  the  coefficients  of  the  incident 
and  scattered  fields  to  the  coefficients  of  the  field  in  Region  1 


/,  =  +2:(Re.O“*K} 

o,=<X{ei,  ,  (Out,Re)ai-  +  (Out,  Out)y9,l- 


where 


2i„.(Re,Re)=  ■  |^[VRev'„(A:or)Re^^/„,(^,r)-  Ret</„(A:or)VRe^«^„,(A:,r)]| 


The  first  argument  (Re  or  Out)  in  Ql„.  signifies  taking  Rev'„  or  \(/„  as  the  wave 
functions  with  argument  k^r  in  the  integral;  the  second  argument  (Re  or  Out) 
signifies  taking  Rey^„  or  y/^  as  the  wave  functions  with  argument  /t,r . 

Next,  we  apply  Huygen’s  principle  to  the  volume  between  o^  and  02.  Expressing 
the  field  inside  02  as 


n 


Using  the  acoustic  boundary  conditions  on  surface  02,  and  evaluating  the 
integrals  yields  two  equations  relating  a\,al,  and  . 

=  ^C™'(Out,Re)a„', 

n' 


23 


Written  in  vector  form,  our  results  for  the  two  layer  scatterer  are 


a  =  /{q'  (Out,Re)a'  +  Q'  (Out,  Out)3' 
f  =  -;{q'  (Re,Re)a'  +  Q'  (Re,  Out)p' } 
a'  =  7Q^(Out,Re)a^ 
p'  =/Q'(Re,Re)a' 


Solving  these  for  the  relation  between  a  and  f 
layered  scatterer 


Jl2  ^  _  Ql  [qI  ||j  Ql 


yields  the  T-matrix  for  a  two- 

(Out,  Out)T^  [q  '  (Out,  Re)] ' 


where 


T'  =-Q'(Re,Re)[Q'(Out,Re)]'' 

is  the  T-matrix  of  a  non-layered  scatterer  bounded  by  surface  ai  and  having  the 
material  properties  of  Region  /. 

Reference  24  give  a  recursive  approach  for  calculating  the  T-matrix  of  an  N  layer 
scatterer.  The  T-matrix  for  a  scatterer  formed  by  removing  the  first  j  layers  of  an 
N  layer  scatterer  layers  ’  is  given  by 


^  _Q7»(Re^Out)T^^' . (Out, Re)]"' 

X  {l  +  Q'"'(Out,Out)T'"' . ^  [q'"'  (Out,Re)]"' } 


which  can  be  constructed  by  working  from  the  innermost  layer  out  to  the  surface. 
E.  Inhomogeneous  Scatterer 

Scattering  from  an  object  with  continuously  varying  properties  can  be  treated 
within  the  T-matrix  formalism,  but  general  results  cannot  be  given.  The  difficulty 
is  that  the  fields  in  the  interior  of  the  object  cannot  be  expanded  in  terms  of 
spherical  wave  functions  since  they  may  not  be  solutions  to  the  wave  equation 
for  the  interior  of  the  scattering  object^®.  A  basis  set  must  be  defined  for  each 
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type  of  inhomogeneity.  Once  the  basis  set  is  known,  the  rest  of  the  T-matrix 
formalism  follows. 

F.  Scattering  from  highly  non-spherical  shapes 

The  development  of  T-matrix  theory  that  we  have  presented  uses  spherical 
coordinates  and  basis  functions  and  works  well  for  scattering  objects  that  are 
nearly  spherical.  There  is,  however,  no  fundamental  restriction  to  treating  objects 
with  more  complicated  shapes.  The  theory  itself  places  minimal  restrictions  on 
the  shape  of  the  scatterer.  The  surface  of  the  scatterer  must  be  continuous  and 
described  by  a  single  valued  function  of  spherical  coordinates.  Scattering  from 
nearly  cylindrical  objects  can  be  treated  using  cylindrical  coordinates  and  basis 
functions.  We  will  not  discuss  cylindrical  coordinates  since  the  theory  follows 
that  for  spherical  coordinates  very  closely.  In  this  section  we  discuss  scattering 
from  highly  non-spherical  objects  such  as  prolate  and  oblate  spheroids  of  high 
aspect  ratio. 

The  results  presented  above  can  be  used  directly,  provided  we  know  how  to 
describe  the  surface  element  in  the  integral  equations  for  the Q  matrix.  For  a 
spheroidal  (prolate  or  oblate)  scatterer,  the  surface  element  is^® 


n(0)cf(T  = 


tg  dr{ey 
r  dO  j 


s\ndddd<l> 


where:  a  and  b  are  the  semi-major  and  semi-minor  axes  of  the  spheroid,  and 


'(«)= 


COS  U  sin  6/ 


a 


For  a  prolate  spheroid  (cigar  shape)  a  >  Z) ;  for  an  oblate  spheroid  (pancake 
shape)  a<b  \  and  for  a  sphere,  a  =  b. 

For  objects  with  large  aspect  ratios  {a»b  or  a  « Z> )  there  are  practical 
problems  in  applying  the  T-matrix  approach.  Larger  numbers  of  terms  in  the 
approximation  of  the  infinite  dimensional  Q  matrix  are  needed  for  convergence 
and  the  matrices  become  ill  conditioned.  Techniques  for  dealing  with  the  ill 
conditioned  matrix  problem  have  been  developed^®  ®®  but  will  not  be  discussed 
here. 
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Reference  29  discusses  a  fundamentally  different  approach  to  applying  the  T- 
matrix  formalism  to  large  aspect  ratio  scatterers.  The  idea  is  to  use  a  coordinate 
system  that  closely  approximates  the  shape  of  the  scatterer  to  reduce  the 
number  of  terms  needed  to  approximate  the  matrices  and  reduce  the  severity  of 
the  ill  conditioned  matrix  problem.  This  approach  is  novel  since  it  is  well  known 
that  the  vector  wave  equation  is  separable  only  in  Cartesian,  cylindrical,  and 
spherical  coordinates.'^  An  analytic  Green’s  tensor,  which  is  used  to  construct 
the  T-matrix,  cannot  be  formulated  in  other  coordinate  systems  since  the  vector 
basis  functions  are  not  orthogonal.  Reference  29  avoids  this  problem  by 
developing  a  Huygen’s  principle  based  on  Betti’s  third  identity^^  that  does  not 
require  the  Green’s  tensor.  They  use  spheroidal  vector  basis  functions  that 
contain  the  spheroidal  counterparts  to  the  spherical  Hankel  functions  and 
spherical  harmonics.  This  approach  is  outlined  below. 

Betti’s  identity  states  that  if  u  and  v  are  solutions  to  the  time  independent  vector 
Helmholtz  equation  in  the  volume,  and  f  (u)  and  f  (v)  are  the  stress  tensors 
associated  with  u  and  v ,  then 

fcfon[t(w)v-t(v)u]=0 


where  the  integral  is  over  the  surface  bounding  a  volume  and  n  is  a  unit  vector 
pointing  out  of  the  volume  and  normal  to  the  surface. 

Huygen’s  principle  is  obtained  by  using  expressions  for  u  and  v  which  describe 
the  scattering  process.  We  bound  the  scattering  object  on  the  outside  by  the 
smallest  spheroid  S+  containing  the  scatterer  and  on  the  inside  by  the  largest 
spheroid  S.  contained  by  the  scatterer.  Using  the  spheroidal  basis  functions 
(given  in  Ref.  29)  we  expand  the  incident  wave  u',  the  scattered  wave  u^  and 
the  wave  inside  S.  u'^  as 


n 

n 


where  are  the  spheroidal  vector  basis  functions.  The  superscript  0  indicates 
that  the  properties  of  the  scatterer  rather  than  the  background  medium  are  used. 
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For  the  surface  bounded  by  S,  the  surface  of  the  scattering  object,  and  S+,  we 
have 

J* dcm  ■  [t(m)-  V  -  T(v)-  u]  =  J don  •  [t(m)-  v  -  f  (v)-  u] 

On  S+  we  have 

n  n 

n  n 

where  t  =  n  •  T  is  the  traction.  Choosing  v  =  and  v  =  Re^?„  and  substituting 
above  results  into  the  two  equations  yields 

-  0„.„a„.  =  J do[i^  («)•  ij/„  -  t(y/„  )•  u  J 

S 

'X  "  I  Vn)-^A 

S 


where  and  t_  are  the  traction  and  displacement  on  outside  surface  of  the 
scatterer.  The  matrix  elements  0„.„  are  related  to  the  orthogonality  and 
normalization  of  the  spheroidal  harmonics. 

We  now  apply  the  boundary  conditions  for  two  solids  in  contact  at  the  surface  of 
the  scatterer:  =  t_  and  =  u_ .  Replacing  the  quantities  on  the  outside  of  the 
scatterer  with  those  on  the  inside  in  the  equations  for  the  incident  and  scattered 
wave  coefficients  results  in 

- 0„.„a„.  =  J do\i_{u)- y/„  -  u_ ] 

S 

s 
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We  next  expand  the  surface  traction  and  displacement  in  term  of  the  spheroidal 
basis  functions  for  the  interior  of  the  scatterer. 


n 

t- 


Substituting  these  expansions  into  the  integral  equations  and  using 

Kn'  =  J 

S 

=  j if£T[t(Re^?„.)- Re^?„  - t(Re^?„)-  Re^^„,] 

5 


we  find 


«'  n' 


Since  the  matrix  O  is  symmetric,  we  arrive  at  the  T-matrix  for  a  spheroidal 
system 


T  =  -ReQ(Q)-' 
Q  =  0  R 
ReQ  =  0'ReR 


The  T-matrix  for  the  spheroidal  system  differs  from  that  for  a  spherical  system  by 
the  inclusion  of  the  O  matrices.  The  non-diagonal  nature  of  O  arises  from  the 
mode  coupling  in  geometries  other  than  spherical. 
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G.  Buried  scatterers 


For  application  to  scattering  from  objects  buried  in  the  seafloor,  we  need  to  know 
how  to  use  the  T-matrix  for  a  scatterer  buried  in  a  half-space.  In  this  section  we 
follow  the  presentation  given  in  Ref.  33. 

3.  Preliminaries 

Figure  4  shows  the  geometry  we  will  consider.  There  are  three  acoustic  media  in 
the  volumes  Vo,  Vi,  and  V2.  The  media  have  densities  pi,  velocities  q,  and  wave 
functions  'Tj,  where  i=1 ,  2,  and  3  indicate  the  media.  The  surface  ai  encloses  the 
scatterer  located  below  the  two-dimensional  surface  ao  separating  volumes  Vo 
and  Vi.  The  surface  ao  has  an  average  depth  of  z  =  0  and  is  confined  between 
two  planes  with  constant  z  coordinates  indicated  by  the  dashed  lines.  The  point 
P  is  the  source  location  in  volume  Vo. 


P. 


Vo,  po,  Co 


Figure  4.  Geometry  for  a  buried  scatterer. 


We  assume  that  the  wave  field  is  harmonic  and  generated  by  a  point  source  at  P 
in  Vo.  The  wave  fields  in  each  of  the  media  satisfy 

=0,  =  1,2,3 

where  =  mlc^  is  the  wavenumber  of  medium  r. 
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We  begin  by  considering  the  integral  equations  for  a  surface  bounded  by  a  finite 
part  of  oo  and  half  of  a  sphere  into  the  lower  half  plane.  Letting  the  radius  of  the 
sphere  go  to  infinity  we  obtain: 


X  (■•)+  jd<7'a,-  [t;  (r')7V(*.|r  -  rl)-  g{*„|r  -  rl)vX(r')] 
0  =  X  (■•)+  fi.  ■  k  {r'yg{k,\r  -  rl)-^fc|r  -  (r')] 


for  r  in  V,, 

(26) 

for  r  outside  Vq 


where  'Fg'  is  the  field  on  ao  as  it  is  approached  from  above.  For  this  problem  we 
will  use  an  expansion  of  the  free  space  Green’s  function  in  terms  of  plane  wave 
traveling  away  from  oo  instead  of  the  spherical  wave  expansion  used  above.  By 
restricting  the  region  of  Vo  that  we  consider  to  those  r  for  which  z>z^  and  the 
region  outside  of  Vo  to  those  r  for  which  r  <  we  can  use  the  following 
expansion  for  the  Green’s  function  in  both  integrals. 


gfck  -  '■1)= j 

0 

kj  =ko//ro  where:  k,,  =A:o(sin^cos^,sin^sin^,cos^),  6  and  ^  are  the  polar  and 
azimuthal  angles  giving  the  spherical  coordinates  of  the  unit  vector  k„  =  kg/^p . 
The  integration  is  over  contour  C+  for  z>  z^  (propagation  upward)  and  over  C.  for 
z<z^  (propagation  downward).  The  contours  are  shown  in  Fig.  5. 


Figure  5.  Integration  contours. 
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4.  Integral  Equations  for  the  Surface  ao 

In  this  section  we  obtain  integral  equations  linking  the  incident  field  and  the  field 
scattered  from  the  plane  surface.  Using  Eq.  26  for  r  in  Vo  and  the  expansion  of 
the  Green’s  function  in  terms  of  plane  waves,  the  scattered  field  in  the  region 
above  can  be  written  as 


2;r 


sin^n 


where 


/(!*.)=  "o  •  VX*(f')]  k,  £  C.  (27) 


We  next  use  Eq.  26  for  r  outside  Vo  to  obtain  an  expression  linking  the  incident 
field  to  the  fields  on  the  surface.  The  incident  field  emanating  from  a  point 
source  located  at  in  Vo  above  is  expressed  first  in  terms  of  spherical  wave 
basis  functions  as 


n 


where  the  subscript  n  refers  to  the  subscript  set /?«/w  used  in  previous  sections. 
For  z  >  Zp  and  z  <  z^  the  spherical  wave  functions  can  be  expanded  in  terms  of 
plane  waves  as"*® 


V.  W  =  ^  S'"  » 

0 

where  the  contour  is  over  C+  for  z  >  0  and  over  C.  for  z  <  0 .  For  z  <  z^  the 
expression  for  the  incident  field  in  terms  of  the  plane  wave  expansion  is 
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We  now  use  Eq.  26  for  r  outside  Vo  to  obtain 


fi.  ■  K*(r')V'(e-'''-'')-e-*"’''V'P;(r')]  K  ^ C.  (28) 


The  only  difference  in  the  integrals  for /(ko)  and  a(ko)  is  that  the  plane  wave 
expansions  refer  to  upward  propagating  waves  (kg  eCJfor/(ko)  and 
downward  propagating  waves  (k^  e  C_ )  for  a(ko). 

5.  Surface  Fields 

We  will  now  use  the  integral  equations  based  on  the  surface  bounding  medium  1 
to  remove  the  surface  fields  from  the  expressions  for /(kg)  and  a(ko).  The 
surfaces  oo  and  oi  enclose  Vi.  Since  there  are  no  sources  in  Vi,  the  integral 
equations  for  this  surface  are 

■P,  (r)  =  n„  ■  [>P,- |r  -  r'|)- |r  -  r'|)v'>P,-  (r ')]  for  r  in  V, 

O'o 

+  |rfcr'n,  ■  ['P;(r')V'g(*,|r-r'|)-g(i,|r-rl)5’'>i';(r'^ 

(29) 

0  =  -  ii,  ■  [>p-(r’)V'^(t,|r  -  r’|)-^(i,  |r  -  rl)7T,-(r')]  for  r  in  V„  or  V, 

o'a 

+  |cf(T'n,  •  ['T,^(rOv'g(^,|r-r'|)-g(*Jr-r'|)v'^;(r')] 


The  boundary  conditions  on  ao  and  oi  will  be  used  to  relate  the  surface  fields  'P* 
to  the  surface  fields  and  'P; .  To  eliminate  'P;  we  use  two  forms  of  Eq.  29, 
one  for  r  in  Vo  and  the  other  for  r  in  V2. 

For  r  in  Vo,  the  T-matrix  based  on  spherical  waves  gives  the  effect  of  the 
scatterer  bounded  by  ai.  The  appropriate  Green’s  function  for  the  integral  over 
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oi  is  given  by  the  expansion  in  spherical  basis  functions  as  used  in  the  previous 
sections, 


g(A,  |r  -  r'l)  =  /I,  ^Re  (A,rJ 


To  enable  a  consistent  definition  of  r^and  r^for  the  integral  over  oi,  we  restrict  r 
to  the  exterior  of  a  sphere  with  a  center  inside  ai  and  enclosing  oi.  Then  for  the 
integral  over  the  surface  of  oi,  we  have  =  r  and  =  r' .  Equation  29  for  r  in 
Vo  is  then 


Jrfo-’  n.  ■  [r;  |r  -  r'|)-  g(k,  |r  -  rl)^X(r')l 


(30) 


Similarly,  by  restricting  r  to  the  interior  of  a  sphere  with  a  center  inside  oi  and 
enclosed  by  ai  we  have  =  r'and  =  r.  We  then  have  for  r  in  V2 

\d<7'  lio  •  [^r(r')V'^(fc,|r-r'|)-g(A:,|r-r'|)v'T',-(r')]  (31) 

0-0 

=  »,2]Rer.fer)jrf<T’n,  •  ['J';(r')7V.(i,r')-n(i,r')^'^,"(r'^ 


We  next  use  the  boundary  conditions  to  eliminate  in  the  above  two 
equations.  Since  we  want  to  describe  the  scattering  object  in  terms  of  its  T- 
matrix,  we  first  expand  the  surface  field  inside  the  scatterer  we  expand  in  regular 
basis  functions  and  obtain 


%  (r) = 5]"?^  fe*") 
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For  the  field  in  volume  Vi,  we  use  the  plane  wave  expansion  containing  both 
waves  propagating  both  upward  and  downward.  This  corresponds  to  both 
integration  contours  and  results  in 


2ff 

'F,-  (r)  j de,  a(k, )  e*'"  sin  ^,  +  j dO,  y0(k, ) e'"-"  sin  9, 


Using  the  acoustic  boundary  conditions  on  ai,  Eq.  30  for  r  in  Vo  becomes 


^d<|)^  a(k,)/(k,,r)sin y^(ki)/(k,,r)sin^, 

0  |_c_  c» 


(32) 


and  Eq.  31  for  r  in  V2  becomes 


In 

^dO^  a(k,)/(k,,r)sin0i  +^d6^  y9(k,)/(k,,r)sin 

0  c_ 


=  Wn  (^1  (Out,  Re)al: 


(2) 


(33) 


We  note  that  Eqs.  32  and  33  differ  in  two  ways.  First,  the  arguments  of  the  Q 
matrices  are  different.  Second,  the  field  points  r  are  different.  For  Eq.  32,  r  is 
in  Vo  while  for  Eq.  33  r  is  in  V2. 

The  Q  matrices  are  defined  by 


0w(Re,Re)=  A:,  jifo-'  n,  •  |Rev'„(A:,r')V'Re<^^„,(^2r')--^[V'Rei/„(^,r')]Re^„,(fe2r')| 
6™'(Out,Re)=  A:,  Jfifo-'  n,  •  |?<'„(A:,r')V'Re^/„,(A:2r')-^[VV„(A:,r')]Re^;^„,(A:2r')| 
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and  /(k,r)  is  given  by 


/(k,r)  =  pa’  a.  ■  |r  -  r'|)-  (v’e"‘'')|r(i,|r  -  rl)] 

0-0 


The  value  of  /(k,r)  depends  on  whether  k  is  on  contour  C+  or  C.  and  on 
whether  r  is  above  or  below  oo.  Evaluating  /(k,r)  yields 


/(k.r)= 


/(k,r)  = 


for  k  e  1 
fork  e  C_J 

fork  G cA 
for  k  G  C_  J 


for  r  above  cx^ 
for  r  below  cXq 


Using  these  values  for  /(k,r)  Eqs.  32  and  33  become 


2;r 

=-/^t^„(A:,r)0„„,(Re,Re)aJ)  forrin  (34) 

0 

2n 

^d<j>,^de^  a(k,)(-e'‘‘''^)sin^,  = „{K^)Q„„{OnX,KQ)a^A  forrin  Vj  (35) 

0  C_ 


We  now  introduce  a  quantity  2(ko,k,)  for  scattering  of  plane  waves  from  the 
surface  oo  that  is  analogous  to  the  Q  matrix  for  scattering  from  a  bounded 
scatterer. 


The  value  of  |Q(ko,k,)  depends  on  the  contour  used  for  k,,  and  k, . 

Using  the  plane  wave  expansion  for  'F'  and0(k(,,k,),  Eqs.  27  and  28  for  /(kj 
anda(ko)  become 
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In: 

/(ko)  =  -'J^^i  a(k,)0(ko,k,)sin^,  +  /?(k,)2(ko,k,)sin^, 


0  C_ 


2n: 

^*(ko)='|^^i  «(k,)g(ko,k,)sin^,  +|c/^,  y?(k,)0(ko,k,)sin 


0  C. 


ko  e  C,  (36) 


koeC  (37) 


6.  T-matrix  for  the  Buried  Scatterer 

Equations  34  and  35  involve  expansions  using  two  different  basis  sets.  The  left 
side  of  these  equations  is  an  expansion  in  plane  waves  while  the  right  sides  are 
expansions  in  spherical  harmonics.  These  equations  are  further  manipulated 
and  used  to  express  the  relationship  between  a(k,)  and  y9(k,)  in  terms  of  the  T- 
matrix  for  the  scattering  object.  First  we  manipulate  Eq.  35  into  a  relationship 
based  on  the  spherical  wave  function  basis  set.  Multiplying  each  side  by  the 
spherical  harmonic  T^(r),  integrating  over  the  unit  sphere,  and  using  the  plane 
wave  expansion  of  the  spherical  wave  functions  on  the  unit  sphere 


2n  n 

jm  (r)  =  j dcp^ dS  (f )sin  & 

0  0 


yields 
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-  (^)  j d<t>,  a(k , )  4  (k ,  )sin  9, 
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27t  Jt 


jd<pjdSY.{r)Y,(r)sm9 


0  0 


0„„.(Out,Re)aJ 


(2) 


The  quantity  in  brackets  above  is  the  orthogonality  integral  for  spherical 
harmonics  and  is  S„  „ .  Using  this  we  have 
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\d(t>,^de,  a(k,)7„(k,)sin^,  =/^2„„,(Out,Re)aP 
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(38) 
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Equation  34  can  be  similarly  manipulated  to  give  a  relationship  between  plane 
wave  amplitudes.  By  considering  r  on  a  plane  of  constant  z  =  Zq  >  we  can 
use  the  plane  wave  expansion  on  contour  C+  for  y/ „{}(/)  and  obtain 


^d(f),^de,  y0(k,)(e''‘'")sin(9, 
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=  j j ^i2„„'(Re,Re)a? 
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This  equation  is  true  provided 

nn 


Equations  38  and  39  can  be  further  manipulated  to  introduce  the  T-matrix  of  the 
buried  scatterer.  First  we  solve  Eq.  38  for  the  expansion  coefficients  of  the 
spherical  wave  basis  functions  by  multiplying  both  sides  by  the  inverse  of  the  Q 
matrix  to  give 


^4;r/" Iq-'  (Out,Re)L  J d^,  j d0,  a(k, )  Y„ (k,  )sin  0, 

"  0  C_ 

=  (Out,  Re)]„„  Q„„.  (Out,  Re)ai?^ 


Using  the  fact  that 


X  b"'  (Out,Re)  =  S„„. 


we  have 
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Substituting  this  into  Eq.  37  and  using  the  definition  of  the  T-matrix  of  the 
scattering  object  r„„(l) 


(l)  =  X!  2-.'  b  ■'  (Out,  Re)L 


yields  the  relation  between  the  upward  propagating  and  downward  propagating 
plane  wave  amplitudes  in  Vi. 


2k 

0  c_ 

Equation  40  is  the  analog  to  the  equation  relating  the  incoming  and  scattered 
spherical  wave  amplitudes,  f  =  Ta ,  that  we  derived  above  for  acoustic  scattering. 
Here,  we  have  a  more  complicated  equation  since  we  are  dealing  with  plane 
wave  amplitudes  and  the  spherical  wave  T-matrix  of  the  scatterer.  The  spherical 
harmonics  appearing  in  Eq.  40  and  the  integration  are  related  to  the 
transformations  between  the  plane  wave  and  spherical  harmonic  representations 
of  the  field. 

In  symbolic  terms  we  can  write  the  three  remaining  equations  describing 
scattering  from  a  buried  object,  Eqs.  36,  37,  and  40,  as 


p  =  Ta 


In  these  equations  we  have  expressly  identified  each  of  the  Q  functions  as  being 
different.  Formally  solving  these  equations  yields 


f  =  -Q«(QiQ/,  +T)Q;!,(l+Q.„TQ;l)''a 


When  the  surface  oo  is  flat,  the  Q  functions  can  be  evaluated  in  terms  of  the 
reflection  coefficient  of  the  surface  to  give 
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€  =  («., +tXi-«„T)-' a 


where:  Roi  is  the  reflection  coefficient  of  a  wave  striking  the  interface  from  above 
and  Rio  is  the  reflection  coefficient  of  a  wave  striking  it  from  below.  The  total  T 
matrix  for  the  buried  object  is  then 


Noting  that  ,  this  equation  can  be  expanded  in  powers  of  the  T  matrix 

of  the  buried  object  to  give  the  series 


Since  =  roiTjo ,  where  Toi  and  Tio  are  the  transmission  coefficients  of  the 

interface,  we  have 


--^01  +^oi^oT  +  T(„7joTi?,oT  +  -*- 


This  expression  has  a  simple  physical  interpretation.  The  T  matrix  of  the  buried 
object  is  made  up  of: 

1 .  A  reflection  from  the  interface  with  no  interaction  with  the  buried  scatterer 

2.  Transmission  through  the  interface,  interaction  with  the  scatterer,  and 
transmission  back  through  the  interface 

3.  Transmission  through  the  interface,  interaction  with  the  scatterer, 
reflection  from  the  underside  of  the  interface,  interaction  with  the  scatterer, 
and  transmission  back  through  the  interface 

4.  ... 

H.  Scattering  in  a  Waveguide 

For  application  to  the  analysis  of  scattering  in  shallow  water,  we  need  the  theory 
of  T-Matrix  scattering  in  a  waveguide.  This  is  developed  in  Refs.  35-37.  Our 
presentation  here  follows  Ref.  37  which  presents  a  compact  theory  for  a 
homogenous  waveguide  that  illustrates  the  basic  physics  without  the 
complications  of  more  general  waveguide  treated  in  Refs.  35  and  36. 

We  assume  a  horizontally  stratified  medium.  The  water  column  is  bounded 
above  by  the  sea  surface  and  below  by  a  layered  seafloor  terminated  by  a  half 
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space.  To  develop  the  basic  concepts,  we  first  consider  the  case  where  both  the 
source  and  scatterer  are  located  in  the  same  layer.  After  that,  we  will  discuss  the 
more  realistic  case  where  the  source  and  scatterer  are  in  different  layers  as 
would  be  the  case  for  a  source  in  the  water  and  a  scatterer  buried  in  the  seafloor. 
Cylindrical  coordinates  (p,  z,  (]))  with  an  origin  inside  the  scatterer  are  used.  The 
vector  p  is  the  position  in  the  x-y  plane.  A  harmonic  time  dependence  (e'"”')  is 
assumed. 

7.  Source  and  Scatterer  in  the  Same  Layer 

The  source  and  scatterer  are  bounded  above  by  the  plane  z  =  a  and  below  by 
the  plane  z  =  -b.  The  total  acoustic  field  is  written  as  the  sum  of  a  source 
component  and  scattered  component  . 


O  =3)  +0 

tot  src  '  ^  scat 


Each  of  the  field  components  is  expanded  in  a  multiple-scattering  series  that 
takes  into  account  the  interaction  with  the  medium  above  and  below  the  layer 
containing  the  source  and  scatterer. 


00  4 


0 


src  /  .  /  ji 

;=0  1=1 


The  4  terms  in  the  sum  over  /  are:  the  direct  path  (/=0),  the  path  reflecting  once 
from  the  upper  boundary  {i=1),  the  path  reflecting  once  from  the  lower  boundary 
(/=3),  and  the  path  reflecting  once  from  each  boundary  (/=4).  The  direct  path 
term  is  given  by  the  free-field  Green’s  function  expanded  in  terms  of  upward  and 
downward  propagating  plane  waves 


where  the  field  point  is  given  by  the  cylindrical  coordinates  (p,  z,  (f))  and  the 
source  position  is  given  by  (ps,  Zs,  <t)s).  The  positive  sign  in  the  exponent  applies 
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when  the  field  point  is  above  the  source  position  {z>Zs)  and  the  negative  sign 
when  the  field  point  is  below  the  source  position  (z>Zs).  The  horizontal 
component  of  the  wave  vector  k  is  q  and  the  vertical  component  is 
h(q)  =  -  q^f^ .  For  simplicity  in  keeping  track  of  the  signs  in  the  exponents, 

we  assume  that  z>Zs>0  in  the  following. 

For  the  /=f  term,  we  have  plane  waves  propagating  downward  from  the  source 
(z=Zs)  to  the  boundary  (z=-b)  followed  by  a  reflection  and  propagation  from  z=-b 
to  the  field  point. 


0 

where  F(q)  is  the  plane  wave  reflection  coefficient  of  the  lower  boundary,  which 
will  depend  on  q  in  the  general  case. 

Similarly  we  have 

0 

00 


where  u(q)  is  the  plane  wave  reflection  coefficient  of  the  lower  boundary  of  the 
layer.  The  i=4  component  is  the  multipath  propagating  downward  at  both  the 
source  and  field  point  and  having  one  interaction  with  each  boundary.  These 
four  terms  contain  the  direct  path  and  three  of  the  four  multipaths  having  one 
interaction  with  each  boundary.  The  fourth  multipath  (propagating  upward  from 
the  source  and  arriving  from  below  at  the  field  point  and  interacting  once  with 
each  boundary)  is  in  the  j=1  part  of  the  sum. 

They=f  terms  in  the  sum  differ  from  the  J=0  terms  by  an  additional  scattering 
from  both  the  upper  and  lower  boundaries.  This  is  accomplished  by  multiplying 
they=0  terms  by  a  factor  UF  exp(2ih(a+b)) .  Note  that  for  the  i=1  term,  the 
multiplication  of  the  direct  path  by  this  factor  produces  the  fourth  multipath  that 
was  missing  from  the  J=0  sum. 
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The  terms  for  any  j  are  obtained  from  the  j=0  terms  by  multiplying  by  the  factor 
[UV  exp{2ih{a  +  Z>))]^ .  The  resulting  sum  for  the  source  field  is 


(41) 
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>=0 


The  geometric  series  can  be  summed  to  give 


00 


The  direct  path  contribution  of  the  scatterer  is  the  T-matrix  of  the  scatterer. 
Using  the  spherical  wave  representation,  this  contribution  is 


pml  p*m'V 


pml,p’m'V  ^  p'm'V 


where:  y^p„,{r)=  h,{kr)Yp„,{^,(^)  is  an  outgoing  spherical  wave,  T  is  the  transition 
matrix  of  the  scatterer  in  spherical  coordinates,  and  Cp„,  are  the  expansion 
coefficients  of  the  total  wave  incident  on  the  scatterer. 

To  use  this  in  the  cylindrical  coordinate  system  of  our  waveguide  we  need  the 
transformation  between  cylindrical  and  spherical  wave  functions.  For  the  z  axis 
of  the  scatterer  parallel  to  the  z  axis  of  the  coordinate  system,  this  transformation 
is  given  by 


f±-T‘ 

kj 


J.igpY"" 


COS/M^, 

sin  m^. 


p  =  even 
p  =  odd 


where  the  plus  sign  is  used  for  z>0  and  the  minus  sign  for  z<0.  The  functions 
convert  the  spherical  functions  into  cylindrical  functions  and  are  given  by 
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where  P”  are  the  associated  Legendre  functions. 

Assuming  that  the  field  point  is  above  the  scatterer  (z>0)  we  use  the  plus  sign  in 
the  transformation  equation  and  obtain  the  following  expression  for  the  scattered 
field  can  be  written  as 


(42) 
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The  terms  in  the  sum  over  j  are  :  the  scattered  field  propagating  directly  from  the 
scatterer  to  the  field  point,  the  scattered  field  reflecting  from  the  lower  interface, 
the  scattered  field  reflecting  from  the  lower  interface,  and  the  scattered  field 
reflecting  once  from  each  interface.  As  above,  the  term  raised  to  the/"  power 
adds  the  additional  multiples  to  the  sum.  The  result  is 
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where  the  vector  A  is  given  by 
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The  next  step  is  to  determine  the  expansion  coefficients  of  the  total  field  incident 
on  the  scatterer.  This  is  done  using  the  “effective  field”  concept  in  which  the 
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incident  field  is  expressed  as  the  source  field  and  the  scattered  field  with  the 
direct  scattering  contribution  removed,  i.  e. 


O 

macent 


Since  none  of  the  three  terms  above  have  singularities  within  the  layer  containing 
the  source  and  scatterer,  they  can  be  expanded  in  terms  of  regular  spherical 
wave  functions.  For  this  expansion,  we  have  z<z^  since  we  are  expanding 
about  the  origin  of  coordinates  located  within  the  scatterer.  Note  that  this 
requires  using  the  opposite  sign  for  the  exponential  in  the  free-field  Green’s 
function  than  we  have  been  using.  For  the  incident  wave  we  have  the  expansion 
used  above, 


^indcenl  (44) 

pml 


For  the  source  term  we  will  convert  Eq.  41  to  spherical  coordinates.  The  first 
step  is  to  separate  out  the  source  and  receiver  parts  of  the  Bessel  function  in  Eq. 
41  using 

00 

-P.|)= 

m=0 


Next,  expand  the  cosine  and  note  that  the  two  terms  are  equivalent  to  a  sum 
over  parity  to  get 


^0  fe|p  -  P.  I)  =  )  1  . 


p  m=0 


[cos  /w(^  )cos  )j 

sin  w(^)sin 


This  expression  is  nearly  that  of  the  regular  cylindrical  basis  functions. 


ReX,.(r) 


cos/w(^)] 
sin  /w(^) 
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The  missing  exponential  is  contained  in  the  multipath  contributions  to  Eq.  41 . 
Finally,  using  the  transformation  of  the  regular  cylindrical  basis  functions  to 
regular  spherical  basis  functions, 


ReX,.(r)=J]£., 

/=0 


yields  the  expression  for  the  source  field 

pml 


Here,  the  vector  A  is  the  same  as  that  defined  above  for  Eq.  42. 

Using  the  same  approach,  the  final  term  in  the  effective  field  equation  can  be 
obtained  as 
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where  R  is  the  rescattering  matrix  given  by 
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Combining  Eqs.  44-46  in  the  equation  for  the  effective  field  and  equating  the 
coefficients  of  the  spherical  basis  functions  gives  us  an  equation  for  the 
expansion  coefficients 


^pml  -  Ppml.p’mT  ^  p’mT 
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p’^mT 
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In  matrix  notation,  this  equation  is  c  =  A+/RTc  which  can  be  solved  for  the 
quantity  Tc  that  appears  in  Eq.  43, 


Tc  =  /(/T-'  +r)''a 

Using  this  in  Eq.  43  yields  the  final  result 

'1>„«=4>„+}a"(;T-'+r)-'a 


This  equation  gives  the  total  acoustic  field  as  the  sum  of  the  source  field  (with 
boundary  interactions)  and  the  scattered  field.  Expanding  the  term  containing 
the  R  and  T  matrices  gives  the  series 


/T  +  TRT + /TRTRT  +  •  • 


The  contribution  of  the  scatterer  to  the  total  field  is  the  direct  scattering  term  with 
corrections  for  the  rescattering  that  occurs  when  the  scattered  field  reflects  from 
the  boundaries  and  is  again  incident  upon  the  scatterer. 

8.  Source  and  Scatterer  in  Different  Layers 

Figure  6  shows  the  geometry  for  source  and  scatterer  in  different  layers.  The 
source  is  in  a  layer  bounded  by  zi  and  Z2.  The  scatterer  is  in  a  layer  bounded  by 
Z3  and  Z4.  The  scatterer  could  be  located  in  the  seafloor  or  in  the  water  column, 
depending  on  the  location  of  Z3. 

For  this  geometry,  the  Green’s  function  contains  the  transmission  coefficient  of 
the  acoustic  field  through  the  layers  between  the  source  and  scatterer.  In  the 
following  ,  we  refer  to  the  layer  containing  the  source  as  layer  1  and  the  layer 
containing  the  scatterer  as  layer  2.  We  will  use  subscripts  1  and  2  to  associate 
the  quantities  h,  k,  U,  and  V  with  the  layer  containing  the  source  or  the  scatterer. 
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Sea  Surface 


Z=Zi - 
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Figure  6.  Geometry  for  source  and  scatterer.in 
different  layers 


The  Green’s  function  for  this  cross-layer  geometry  is 


00 


where:  /,  (zj  is  the  total  acoustic  field  in  layer  2  incident  upon  the  interface  at  Z2, 
f^{z)  is  the  total  field  at  the  point  z  in  layer  2,  and  l/Vis  the  total  plane  wave 
transmission  coefficient  through  the  layers  between  layer  1  and  layer  2.  After 
carrying  out  the  multipath  expansion  of  the  fields  in  layers  1  and  2  as  in  the 
previous  section,  we  find 


\.  -.V 

/itJ  ^ 

,  rr  -.r  2ih.(2.  2A 


Inserting  these  into  the  Green’s  function  yields  the  field  produced  by  the  source 
at  a  field  point  in  layer  2, 


4,  =  _L ]Aj  (,|p - p  IV-V. 
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The  equation  for  A(rJ  for  the  cross-layer  geometry  is 
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21^2(23+24 ) 


The  rescattering  matrix  for  the  cross-layer  geometry  is  obtained  from  that  above 
by  substituting  23  for  a,  Z4  for  -b,  k2  for  k,  U2  for  U,  and  V2  for  V. 

9.  Normal  Mode  Formulation 

Moving  the  integration  contour  into  the  upper  right  quadrant  recovers  the  normal 
mode  formulation  for  the  within-layer  and  cross-layer  problems.  The  poles  of  the 
integrand  occur  at  the  zeros  of  the  denominators.  For  the  within-layer  case  the 
poles  are  defined  by  the  equation  (l  -  0 .  For  the  cross-layer  case 

the  defining  equation  0. 

I.  Summary 

The  T-matrix  approach  is  fully  developed  and  contains  all  the  tools  needed  for 
studying  scattering  from  geological  features  buried  in  the  seafloor.  The  theory  for 
both  acoustic  and  elastic  wave  scattering  exists  as  well  as  extensions  to  multiple 
scattering  and  scattering  in  a  waveguide. 

IV.  P  S  COUPLING  DUE  TO  GRADIENTS 

Independently  propagating  shear  and  compressional  waves  can  exist  only  when 
the  acoustic  properties  of  a  material  do  not  depend  on  position.  When  gradients 
exist,  the  two  types  of  waves  couple  within  the  material  as  well  as  at  interfaces 
between  two  different  materials.  Expanding  in  the  inverse  of  the  frequency, 
Richards^^  developed  an  asymptotic  theory  for  the  coupling  of  shear  and 
compressional  wave  potentials  and  demonstrated  decoupling  at  high 
frequencies.  The  unique  and  useful  aspect  of  this  theory  is  the  use  of  second 
order,  rather  than  fourth  order,  differential  equations  for  the  coupled  potentials. 
Vidmar”*^  extended  this  theory  to  a  form  valid  for  all  frequencies  that  became  the 
basis  for  estimates^  that  gradient-driven  coupling  is  important  only  below  about 
10  Hz  for  deep  sea  sediments.  Numerical  studies^  verified  this  estimate.  Other 
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research^  showed  that  gradient  driven  coupling  generates  shear  waves  when 
normal  coupling  at  an  interface  does  not  take  place,  i.  e.,  when  the  shear  speed 
is  zero  at  the  interface. 

Recent  borehole  tomography  measurements^^  in  shallow  water  indicate  that 
volume  inhomogeneities  have  gradients  that  are  significantly  larger  (10  times,  or 
more)  than  those  found  in  deep  ocean  sediments.  These  enormously  high 
gradients  raise  the  possibility  that  gradient-driven  PS  conversion  could  be  an 
important  acoustical  process  at  frequencies  up  to  1  kHz.  A  recent  theory^  of 
scattering  from  these  inhomogeneities  has  been  carried  out  but  does  not  include 
shear  wave  coupling. 

The  remainder  of  this  section  reviews  the  theory  of  PS  coupling  as  developed  in 
Refs.  42  and  43. 

A.  Theory  for  Homogeneous  Media 

The  description  of  motion  in  terms  of  potentials  for  a  homogeneous  layer 
illustrates  the  fundamental  approach  used  for  heterogeneous  media.  We  begin 
with  the  equation  of  motion  for  the  vector  displacement  u  for  a  material  with 
density  p  and  Lame  parameters  X,  and  p 

= -(A  +  2/^)VV-u  +  //Vx  Vxu  (47) 


and  the  decomposition  of  the  displacement  into  compressional  and  shear 
potentials  as 


u^VO  +  Vx'i'  (48) 


The  usual  procedure  is  to  substitute  Eq.  48  into  Eq.  47  (assuming  that  V  •  'F  =  0 ) 
yielding  a  third  order  differential  equation.  Further  manipulation  results  in  two 
decoupled  fourth  order  differential  equations  for  the  compressional  and  shear 
wave  potentials. 


V' 
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=  0 


49 


where  we  have  used  the  definitions  of  the  shear  speed  =  nj p  and  the 

compressional  speed  c\  ={X  +  2p)lp. 

Upon  assuming  that  the  quantities  in  brackets  are  zero,  we  have  the  usual 
second  order  differential  equations  for  the  potentials.  When  applied  to 
heterogeneous  media,  this  approach  yields  coupled  fourth  order  differential 
equations.'^ 

Another  approach  is  to  note  that  Eq.  47  already  has  the  form  of  Eq.  48  and  to 
then  define  the  potentials  directly  from  Eq.  47  as 


^  T7 

— 5-V  u 
vr'^  V 

- ;-VxU 

c: 


It  follows  from  Eq.  47  that 

v^o+^a)  =  o 

■7  -  — 

v-^  =  o 

This  approach  avoids  using  fourth  order  differential  equations  and,  when 
generalized  to  heterogeneous  media,  it  produces  coupled  second  order 
differential  equations. 

B.  Theory  for  Heterogeneous  Media 

We  begin  with  the  equation  of  motion  for  a  heterogeneous  solid  in  Cartesian 
coordinates^.  For  simplicity,  we  consider  wave  motion  in  the  x-z  plane  in  a 
material  with  gradients  in  the  z  direction  only.  The  time  dependence  is  e'"".  The 
dependence  on  x  is  given  by  e'*  which  follows  from  spatial  invariance  along  the  x 
axis. 
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1 .  Weakly  Coupled  Potentials 
The  resulting  equation  of  motion  for  a  harmonic  wave  is 


pco^vi  =  -(/l  +  2//)VVu  +  //Vx  Vxu-lzVu  -  //' 


r 


zxV  xu  +  2 


V 


dz  ^ 


(49) 


where  z  is  a  unit  vector  in  the  direction  of  the  positive  z  axis  and  the  prime 
indicates  a  derivative  with  respect  to  z. 

We  next  assume  that  the  displacement  can  be  written  in  terms  of  potentials  with 
some  additional  functions  that  will  allow  us  to  recover  Eq.  49. 

u  =  — va)+— Vxi'-s  (50) 

f\  fi 


For  motion  in  the  x-z  plane,  the  vector  potential  is  given  by  ’i'(x,z)  =  ^^(x,^) 
where  y  is  a  unit  vector  along  the  positive  y  axis.  As  in  the  homogeneous  case, 
we  define  the  potentials  in  terms  of  the  displacement  vector,  but  with  some 
additional  (as  yet  unknown)  functions  in  the  definition. 

®  =  aV(/,ii)  (51) 

'i'=>9Vx(/,u) 


The  functions  /, ,  /j ,  a,  and  p  depend  on  z  alone.  They  will  allow  Eq.  50  to 
duplicate  the  first  and  second  derivative  terms  in  Eq.  49.  The  vector  s  contains 
the  remaining  terms  in  Eq.  47  that  are  linear  in  the  components  of  u.  We 
proceed  by  calculating  the  quantity 

u’  =— VO+— Vx^ 

/. 


and  comparing  the  resulting  expression  with  Eq.  49. 
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The  two  terms  of  u*  are 


ivO  =  V(V-  u)+z(a /j' V  •  u  +  -^f z  x  (V  x  u)+— 1  +  feiz(z  •  u)  and 

/i  /i  /i  ^  /i 

J_  V  X  T  ^  V  X  (V  X  u)+^^^z  X  (V  X  u)+ 


/: 


A 


/2 


’z(V.u)-^' 
^  ^  dz 


+  .^^£x(zxu) 
/2  ^  ^ 


Comparing  the  terms  containing  the  second  derivative  of  u  in  u*  with  Eq.  49 
provides  the  following  two  equations. 


/2  _  {^  +  2p) 

fx  po)^ 


h  PO)^ 


(52) 


The  first  derivative  terms  provide  the  following  three  equations. 


Defining  g,  =  ////,  and  =  f^jf^ ,  solving  for  a  and  p  in  Eq.  52,  and  substituting 
into  the  first  two  equations  of  Eq.  53  yields  the  following. 


Sx  = 

gi  = 


A  +  'IjLl 


^P' 


X  +  ju 
P 


X  +  /j. 


\P 

V 

V,  P 


2p' 

X  +  2y£/ 


(53) 
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We  note  that  g,  +  gj  =  p7 P  that  gi  and  g2  depend  only  on  the  properties  of 

the  material.  Substituting  a,  p,  gi,  and  g2  into  the  third  equation  of  Eq.  53 
produces  an  identity.  Once  we  know  g,-  in  terms  of  the  material  parameters,  it  is 
easy  to  calculate  f,.  From  the  definition  of  g,  it  is  the  depth  derivative  of  ln(/]). 
Integrating  over  depth  from  the  top  of  the  sediment,  and  assuming  a  value  of  one 
at  the  sediment  surface,  we  have 


We  next  recover  Eq.  49  by  defining  s  as  containing  the  terms  in  u*  which  are 
linear  in  u. 


s  = 


(55) 


Manipulating  the  coefficients  in  Eq.  55  and  substituting  into  Eq.  50  gives  us 


u  = 


±  VO) + -^  V  X 'F + ^ 

/,  A  "  o'  " 


(56) 


where 


=  -(p'gi  +  Pg'i  +  Mg^g^ )  (57) 

P 

^2  =  +  ^P)  g2  +  (^  +  2/^)g2  +  '^P)g2g2 


The  differential  equation  for  O  is  obtained  by  multiplying  Eq.  56  by  fz,  taking  the 
divergence  of  the  result,  and  multiplying  by  a  to  produce  O  as  defined  in  Eq.  51 . 
Multiplying  Eq.  56  by  fi,  taking  the  curl,  and  multiplying  by  p  produces  the 
differential  equation  for  'P.  The  results  are 
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v’o+te,  -?,)4>'+Af^^ji> =/[!);  +^,(A  -A)k  +/,(A  -mK  (58) 

V’T  +  (g,  -  g,  )>!"  +  =  f,  [-  D\+g,  (D,  -  A  )]»,  +  /,  (D,  -  D,  )< 

\,  M  J 


In  Eq.  58  we  have  defined  two  auxiliary  functions  Di  and  D2  as 


A=1 

A=1 


(59) 


Equation  58  shows  the  high  frequency  decoupling  of  the  potentials  <1)  and  'P. 
From  Eq.  59,  the  terms  involving  D^  -D^  and  the  derivatives  of  Di  and  D2  are 
proportional  to  ijco^ .  Using  Eq.  56  to  substitute  for  the  components  of  the 
displacement  in  Eq.  58  produces  two  terms  proportional  to  l/o)^  involving  both 
potentials  (hence,  the  coupling  between  them)  and  two  additional  terms 
proportional  to  ijco^  involving  the  components  of  the  displacement.  Continuing 
the  process  generates  a  series  involving  l/o^  to  higher  powers  multiplying 
increasingly  high  derivatives  of  the  potentials.  To  show  decoupling,  one  argues 
that  for  sufficiently  high  frequencies  the  right  side  of  Eq.  58  is  negligible 
compared  to  the  left  side.  Hence,  the  potentials  are  essentially  decoupled  at  high 
frequencies.  While  this  argument  is  true,  it  is  not  clear  just  how  high  in  frequency 
one  must  go  to  achieve  decoupling  in  view  of  the  increasingly  higher  order 
derivatives  of  the  potentials  appearing  on  the  right  side  of  Eq.  58.  We  address 
this  difficulty  below  and  derive  a  quantitative  definition  of  the  decoupling 
frequency.  Decoupling  also  occurs  for  special  materials  for  which  D[  =  Dj  =  0 
and  Z),  =  £>2 . 

2.  Exact  Coupled  Potentials 

We  follow  Ref.  43  to  derive  second  order  differential  equations  for  O  and  'P  that 
are  exact,  i.  e.,  they  do  not  contain  the  components  of  the  displacement.  We  start 
out  by  writing  Eq.  56  in  component  form  and  combine  terms  to  obtain 


D,m 


X 


D^u 


z 


7(V4>).+^(Vx'!'), 

.J\  Jl 

,J\  Jl 


(60) 
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We  now  write  Eq.  60  in  vector  form  by  defining  a  coupling  tensor  C  and  its 
inverse  as 


C  = 


C  '  = 


f  Z),xx 


V 

^  1 


0 


A 


-XX 


0  ^ 
D^iij 

0  1 


(60) 


Substituting  into  Eq.  58  gives  an  equation  for  the  displacement  in  terms  of  the 
potentials  without  extra  displacement  components  on  the  right  side. 


u  =  C-' 


( \  \  A 

— VO  +— Vx^ 
^f\  fl  J 


(62) 


The  next  step  is  to  use  Eq.  62  and  Eq.  51  to  again  derive  the  differential 
equations  for  the  potentials.  As  in  the  previous  section,  we  multiply  Eq.  62  by  f2, 
take  the  divergence,  and  multiply  by  a  to  produce  the  differential  equation  for  O. 
Multiplying  Eq.  62  by  fi,  taking  the  curl,  and  multiplying  by  p  produces  the 
differential  equation  for  'F.  Recalling  that  the  x  dependence  of  the  potentials  is 
through  e'*" ,  the  result  is 
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V3 

(  2  > 
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/2 
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La 

The  terms  on  the  right  side  of  Eq.  63  provide  the  coupling  between  the  potentials 
without  reference  to  components  of  the  displacement. 
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3.  Displacement  and  Stress  for  Exact  Coupled  Potentials 

We  need  to  know  how  to  express  the  components  of  the  displacement  vector 
and  the  stress  tensor  in  order  to  carry  out  numerical  computations.  The  results 
for  the  displacement  follow  directly  from  Eq.  62. 
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u,  =• 
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iS 
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<D-— 

\f\  /a  y 
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The  expressions  for  the  stress  tensor  are  also  straightforward,  but  tedious,  to 
derive.  The  results  are 


A  =  -p\ 
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(O'^- 
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A  =  P 
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2  ) 


\  O  g2c\  O'  iclW  T' 

J/  "  Sx  ”  A  £1,  TT 

I.*  ■  g,cliS  9  .  .  Ic^iS 
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4.  Decoupling  of  Exact  Potentials 

Decoupling  of  <I)  and  T  occurs  under  two  conditions.  First,  there  is  decoupling  at 
high  frequencies.  Because  coupling  terms  are  proportional  to  the  derivatives  of 
Di  and  D2  or  to  A,  -D^  decoupling  occurs  at  high  frequencies  in  agreement  with 
Ref.  42.  The  second  kind  of  decoupling  occurs  for  particular  media  for  which 
A'  =  A'  =  0  and  D^=D^.  This  is  the  case  for  the  material  found  in  Ref.  50  for 
which  //  =  //,  (1 + Z>z)^ ,  A  =  A,  (1 + bzf ,  A,  =  //, ,  and  p  =  p^.  Direct  substitution  into 
Eqs.  57  and  59  yields  cof  =col  =  -6p^b^  and  verifies  decoupling  at  all 
frequencies.  Note  that  in  evaluating  Eq.  57  we  do  not  need  to  know  fi  and  f2  but 
only  gi  and  g2,  which  are  defined  in  terms  of  material  parameters. 

5.  Resonance  Frequencies 

Equation  63  has  an  interesting  feature.  There  is  a  potential  resonance  in  the 
system  at  the  frequencies  ©1  and  ©2  for  which  the  denominators  in  several  terms 
of  Eq.  63  can  be  zero.  The  significance  of  this  resonance  has  not  been  fully 
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explored.  Here,  we  examine  the  resonant  frequencies  as  a  function  of  the 
gradients  of  seafloor  parameters. 

The  first  step  is  to  develop  expressions  for  the  resonant  frequencies  based  on 
density,  shear  velocity,  compressional  velocity.  Equation  57  gives  the  resonant 
frequencies  in  terms  of  gi,  g2,  the  density  and  the  depth  derivative  of  the  Lame 
parameters.  Our  goal  here  is  to  develop  an  expression  that  involves  more 
readily  obtainable  quantities  such  as  density,  shear  velocity  and  compressional 
velocity.  We  start  by  introducing  the  shear  and  compressional  velocities  into  Eq. 
57  and  eliminating  the  derivatives  of  the  Lame  parameters  using  the  following 
relations  obtained  from  the  derivative  of  the  definition  of  the  shear  and 
compressional  velocities 

M  P 

{X  +  2^)  _p'  ^ 

1  +  2^  p 


The  resulting  resonant  frequencies  are 
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gi,  and  g2  and  their  derivatives  are  given  by 
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a)  Resonance  Frequencies  for  Typical  Sediments 

In  this  section  we  estimate  the  resonance  frequencies  for  typical  sediments.  First, 
we  will  use  typical  values  for  gradients  to  identify  the  driving  terms  in  the 
expressions  for  gi,  g2,  and  their  derivatives.  Typical  parameters  for  near  surface 
sediments®^  are  given  in  Table  2. 


Table  2.  Typical  values  of  surficial  sediment  parameters. 


Using  the  fact  that  clfc\  » IxlO'^  for  clay  and  assuming  that  second  derivatives 
are  negligible,  the  expressions  for  gi,  gz,  and  their  derivatives  are  approximately 
given  by 
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We  next  identify  the  dominant  terms  involving  the  parameter  gradients.  For  deep 
ocean  sediments  (clay),  typical  gradients  are:  p'  =  1.2x10“^  g/cc/m  ,  =  1.5  /s , 

and  c'  =  4.0  /s .  The  values  of  the  ratios  of  the  parameters  to  their  derivatives  are 
approximately  p'/p »lxlO"Vm,  c'JCj,  «lxlO'Vm,  and  c'Jc^  »lxl0''/m.  The 
relative  gradients  of  density  and  compressional  velocity  are  comparable  while  the 
relative  gradient  of  shear  speed  is  about  two  orders  of  magnitude  larger. 

Keeping  only  the  most  significant  terms  in  the  expressions  for  the  resonance 
frequencies,  we  obtain 


«  2cl  and  col  =  4c] 
P 


Substituting  numerical  values  we  estimate  the  resonant  frequencies  to  be 


/  «-^«.07Hz  and  ^«1.6Hz 
'  2;r  iTT 


Thus,  gradient  driven  coupling  of  compressional  and  shear  waves  is  negligible 
for  deep  sea  sediments  for  frequencies  above  about  10  Hz.  Since  parameter 
values  for  shallow  water  sediments  do  not  change  by  orders  of  magnitude,  the 
same  conclusion  is  true  for  shallow  water  sediments. 

b)  Resonance  Frequencies  for  Volume  Inhomogeneities. 

As  discussed  at  the  beginning  of  this  section,  there  is  evidence"*^  that  volume 
inhomogeneities  in  shallow  water  sediments  have  parameter  gradients  that  are 
significantly  larger  (10  times,  or  more)  than  the  large  scale  gradients.  In  this 
section,  we  present  estimates  of  the  resonance  frequencies  associated  with 
these  gradients.  To  estimate  gradients  we  assumed  a  transition  from  a  clay  to  a 
sand  sediment  over  distances  from  1  cm  to  20  m. 

Figure  7  shows  that  the  resonance  frequency  can  be  quite  large  for  a  thin 
transition  layer.  Frequencies  above  1  kHz  occur  for  transition  layers  of  1-5  cm, 
not  an  unreasonable  thickness  for  volume  inhomogeneities.  For  transition  layers 
about  10m  thick,  the  resonance  frequencies  are  on  the  order  of  1  Hz,  the  value 
estimated  for  deep  water  sediment  types. 
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Figure  7.  Dependence  of  resonance  frequency  on  transition  layer  thickness. 
C.  Summary 

We  reviewed  the  theoretical  treatment  of  gradient  driven  coupling  of 
compressional  and  shear  waves.  The  theory  of  a  homogeneous  layer  showed 
the  basic  approach  for  developing  second  order,  rather  than  fourth  order, 
differential  equations  for  the  coupled  potentials.  The  basic  concept  introduced 
was  that  of  defining  the  potentials  in  terms  of  the  displacements.  Applying  this 
concept  to  a  depth  dependent  layer  yielded  two  differential  equations  for  weakly 
coupled  potentials,  for  which  the  coupling  is  proportional  to  the  inverse  of  the 
square  of  the  frequency.  The  potentials  reduce  to  the  ordinary  shear  and 
compressional  potentials  for  a  homogeneous  media.  The  two  potentials  are 
coupled  by  terms  in  the  displacements.  A  refinement  of  the  theory,  using  tensor 
scale  factors,  produced  exact  coupled  differential  equations  with  the  coupling 
terms  related  to  the  potentials  and  their  derivatives. 

The  differential  equations  for  the  exact  coupled  potentials  demonstrate 
decoupling  at  high  frequencies  and  at  all  frequencies  for  media  with  specific 
properties.  Decoupling  at  high  frequencies  occurs  when  the  frequency  is  large 
compared  to  two  frequencies  whose  values  depend  on  the  material  parameters 
and  their  gradients.  These  frequencies  are  called  “resonance  frequencies”  since 
they  appear  in  the  denominators  of  terms  in  the  differential  equations.  Should 
the  frequency  be  equal  to  one  of  the  resonance  frequencies,  the  denominators 
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vanish.  The  significance  of  these  resonance  effects  has  not  been  examined.  For 
typical  deep  sea  sediments  the  resonance  frequencies  are  around  1  Hz.  They 
can  be  on  the  order  of  1 000  Hz  for  the  large  gradients  thought  to  occur  for 
volume  inhomogeneities  in  shallow  water  sediments. 

V.  RECOMMENDATIONS 

Our  investigation  of  T-Matrix  theory  and  the  theory  of  gradient  driven  coupling  of 
shear  and  compressional  waves  finds  that  both  hold  great  promise  for  developing 
an  understanding  of  scattering  from  the  subbottom  in  shallow  water  areas. 

Below  we  give  our  recommendations  for  further  work. 

A.  T-Matrix 

The  T-matrix  approach  to  scattering  is  well  developed  and  contains  almost  all  the 
tools  needed  to  address  some  of  the  issues  dealing  with  the  time  spread  induced 
by  scattering  from  the  seafloor.  One  such  issue  is  the  role  of  scattering  from 
sedimentary  layers  with  finite  lateral  extent.  Scattering  from  the  edges  of  layers 
is  one  possible  mechanism  for  producing  backscatter  from  regions  with  a  smooth 
sediment  surface.  Another  is  the  problem  of  predicting  the  attenuation  parameter 
of  sedimentary  materials  for  Navy  databases.  At  present  the  attenuation  is  a 
“total  energy”  attenuation  which  includes  effects  of  scattering  as  well  as  intrinsic 
attenuation.  Understanding  these  issues  would  provide  the  basis  for 
improvements  to  the  geoacoustic  databases  used  by  the  Navy  to  predict  the 
performance  of  sonar  systems. 

The  approach  we  recommend  would  be  to  first  develop  numerical  models  of  the 
T-matrix  of  several  canonical  objects  that  would  serve  as  entries  in  a  catalog  of 
scatterers.  The  objects  would  include:  (1)  nearly  spherical  layered  objects  that 
would  approximate  a  density  or  velocity  inhomogeneity,  (2)  scatterers 
approximating  the  shape  of  gas  trapped  in  sediments  (spherical  to  disk  shaped), 
(3)  pancake  shaped  (oblate  spheroid),  fluid,  object  with  a  thickness  and  lateral 
extent  typical  of  sedimentary  layering,  and  (4)  objects  approximating  the  shape  of 
other  buried  objects  such  as  rocks  and  shells.  The  numerical  models  must  be 
broadband  so  the  impulse  response  of  the  scatterer  can  be  obtained  as  an 
inverse  Fourier  transform.  Once  the  T-matrix  is  available,  the  methods 
summarized  above  can  be  applied  to  study  scattering  from  single  objects  buried 
in  a  background  sediment.  Next,  scattering  from  several  objects  in  a  typical 
geometry  would  be  studied  to  determine  the  importance  of  multiple  scattering, 
edge  effects,  and  orientation  relative  to  the  horizontal.  Scattering  from  several 
realizations  of  the  layering  would  be  carried  out  to  understand  the  variability  of 
scattering  and  time  spread.  Finally,  scattering  in  a  waveguide  would  be 
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examined  using  simplified  models  based  on  the  earlier  work  with  the  goal  of 
determining  the  importance  of  scattering  in  relation  to  attenuation  and  other 
parameters. 

B.  Gradient  Driven  Coupling 

The  theory  for  gradient  driven  coupling  of  compressional  and  shear  waves  is  well 
developed  with  two  exceptions.  First,  a  stochastic  theory  does  not  exist.  Such  a 
theory  is  needed  to  explore  the  role  of  gradient  driven  coupling  for  volume 
inhomogeneities.  Second,  the  significance  of  the  resonance  frequencies 
described  above  has  not  been  explored.  If  the  resonance  leads  to  instability  in 
the  system,  it  might  be  related  to  the  onset  of  seismic  activity.  For  the  high 
gradients  estimated  for  volume  inhomogeneities  in  shallow  water,  the  resonance 
could  be  an  important  loss  process  converting  compressional  energy  to  shear 
energy  with  greater  efficiency  than  coupling  at  an  interface  between  two  layers. 
We  recommend  further  research  to  address  these  issues. 

A  fully  developed  theory  of  gradient  driven  coupling  would  provide  the  basis  for 
determining  the  importance  of  coupling  to  shear  waves  through  scattering  from 
volume  inhomogeneities.  This  is  an  important  issue  to  resolve  since  it  has 
implications  for  the  content  of  Navy  geoacoustic  databases  and  the  complexity  of 
Navy  performance  prediction  models.  If  shear  wave  coupling  is  important,  new 
propagation  models  will  be  needed  that  will  typically  require  an  order  of 
magnitude  more  computational  resources  than  current  models.  If  shear  coupling 
is  not  important,  simplified  models  and  databases  can  be  used. 
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